---
title: "Replication document"
author: "Francesco Bailo"
date: "30/09/2019"
output: 
  pdf_document:
    toc: true
    toc_depth: 3
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, 
                      message = FALSE, warning = FALSE, fig.path = "figure/", dev = 'pdf',
                      cache = TRUE)

knitr::opts_chunk$set(fig.process = function(x) {
  x2 = sub('-\\d+([.][a-z]+)$', '\\1', x)
  if (file.rename(x, x2)) x2 else x
})

# Colour palette
two_colours <- c("#d7191c","#abdda4")
```

```{r stargazer-set}
sg_param_appendix <- 
  list(omit.stat = "f",
       report = "vc*",
       model.numbers = TRUE,
       model.names = TRUE,
       float = FALSE,
       style = 'ajps',
       star.cutoffs = c(0.05, 0.01, 0.001))
sg_param_intext <- 
  list(omit.stat = c('adj.rsq', 'f', 'ser'),
       report = "vc*",
       model.numbers = FALSE,
       model.names = FALSE,
       float = FALSE,
       style = 'ajps',
       star.cutoffs = c(0.05, 0.01, 0.001))
```

# How to cite the replication materials

```{BibTex}
@data{DVN/OL3EVW_2019,
author = {Bailo, Francesco},
publisher = {Harvard Dataverse},
title = "{Replication Data for: Connecting through the void: Distress, participation and the electoral success of the Five Star Movement}",
year = {2019},
version = {V1},
doi = {10.7910/DVN/OL3EVW},
url = {https://doi.org/10.7910/DVN/OL3EVW}
}
```

# R Packages

```{r}
library(ggplot2)
library(gridExtra)
library(scales)
library(dplyr)
library(knitr)
library(zoo)
library(cowplot)
library(sf)
library(ggsn)
library(eurostat)
library(ggExtra)
library(scales)
library(stringr)
library(stargazer)
library(scales)
library(reshape2)
library(RColorBrewer)
library(readxl)
library(rgeos)
library(sp)
library(plm)
library(viridis)
library(MatchIt)
library(haven)
library(lm.beta)
library(kableExtra)
```

# Variables: computational abbreviations and labels

```{r, echo = FALSE, results = 'asis'}

var_dict_labels = 
  c("pop_density" = "Population density",
    "log(pop_density)" = 'Population density (log)',
    "dist_10k_ppl" = "Distance to first +10,000 pple. city",
    "log(dist_10k_ppl)" = "Distance to first +10,000 pple. city (log)",
    "north" = "North (bin)",
    'unemployment' = "Unemployment (pct)",
    'employment' = "Employment (pct)",
    'partwf_perc' = "Workforce participation (pct)",
    'foreignpop_perc' = "Foreign population (pct)",
    'foreignpop_africa_perc' = "Foreign African population (pct)",
    'over65_perc' = "Population over 65 (pct)",
    'housewive_perc' = "Housewive among population (pct)",
    'degree_perc' = "With university degree (pct)",
    'degree' = "With university degree (bin)",
    "meetup_treated_2018" = 
      "+2 weeks with +1 event (2018)",
    "meetup_treated_2013" = 
      "+2 weeks with +1 event (2013)",
    'meetup2013_5km_buff_freq_7days' = 
      "Weeks with +1 event before 2013 election",
    'meetup2018_5km_buff_freq_7days' = 
      "Weeks with +1 event before 2018 election",
    'treated' = "+2 weeks with event before election (bin)",
    'M5S_2013' = "M5S vote 2013 election (pct)",
    'M5S_2018' = "M5S vote 2018 election (pct)",
    'M5S_2013_neigh_avg' = 
      "M5S vote in neigh. 2013 election (pct)",
    'M5S_2018_neigh_avg' = 
      "M5S vote in neigh. 2018 election (pct)",
    'M5S_diff' = 
      "M5S vote diff. between 2018 and 2013 (pct points)",
    'nat_pol_trust' = "Trust in national institutions",
    'loc_pol_trust' = "Trust in local institutions",
    'pol_trust' = "Overall trust in institutions",
    'income' = "Income (procapite, 10K €)",
    'donor_perc' = "Blood donors (pct)",
    'vol_perc' = "Volunteers (pct)",
    'ANPIV_perc' = "Volunteers (pct)",
    'ref16_turnout' = "Turnout 2016 referendum (pct)",
    'male' = "Male (bin)",
    'edu' = "Education level",
    'age' = 'Age',
    'unemployed' = "Unemployed (bin)",   
    'volon' = "Volunteer of non profit org. (bin)",
    'atgra' = "Volunteer of profit org. (bin)",
    'trust' = "General trust in people (bin)",
    'trust2' = "Turst in security services",
    'trust3' = "Trust in strangers",
    "type_comunePeriphery of metropolitan area" = 
      "Comune type: Periphery of metropolitan area",
    "type_comuneComune up to 2000" = 
      "Comune type: up to 2000",
    "type_comuneComune between 2001 and 10,000" = 
      "Comune type: between 2001 and 10,000",
    "type_comuneComune between 10.001 and 50,000" = 
      "Comune type: between 10.001 and 50,000",
    "type_comuneComune over 50,000" = 
      "Comune type: over 50,000",
    "unemployment:north" = "Unemployment X North",
    "unemployedTRUE:north" = "Unemployed X North",
    "meetup_treated_2013TRUE:nat_pol_trust" = "+2 w. +1 e. X trust nat. inst.",
    "meetup_treated_2013TRUE:eco_stress" = "+2 w. +1 e. X HH eco. distress",
    "meetup_treated_2013TRUE:eco_stress_inv" = "+2 w. +1 e. X HH eco. comfort",
    
    "meetup_treated_2013TRUE:internet_mediated_part" = "+2 w. +1 e. X Internet-med. part.",
    "nat_pol_trust:internet_mediated_part" = "Trust nat. inst. X Internet-med. part.",
    "meetup_treated_2013TRUE:nat_pol_trust:internet_mediated_part"= "+2 w. +1 e. X trust nat. inst. X Internet-med. part.",
    
    "internet_access" = "Internet access (pct of HHs)",
    "eco_stress" = "Household perceived economic distress",
    "eco_stress_inv" = "Household perceived economic comfort",
    "eco_perception" = "National economy perception",
    "poltalk_freq" = "Frequency political talks",
    "online_pol_discussion" = "Frequency online political talks",
    "info_internet" = "Internet primary info. source",
    "info_internet_news" = "News websites primary info. source (bin.)",
    "info_social_media" = "Soc. media primary info. source (bin.)",
    "info_personal" = "Family/Friends primary info. source (bin.)",
    "info_written" = "Newspapers/Magazines primary info. source (bin.)",
    "info_tv" = "TV primary info. source (bin.)",
    "internet_mediated_part" = "Internet-mediated participation (bin.)",
    "interpers_trust" = "Interpersonal trust",
    "member" = "Member of an association (bin.)",
    "pol_interest"  = "Political interest",
    "use_facebook" = "Facebook use (bin.)",
    "female" = "Female",
    "age" = "Age",
    "AmpQ" = "Size of comune",
    "income_norm" = "Income (procapite/family, norm.)",
    "clsness_M5S" = "Closeness to M5S identity",
    "pol_part" = "Political participation",
    "nat_pol_trust:clsness_M5S" = "Trust nat. inst. X closeness to M5S",
    "eco_stress:clsness_M5S" = "HH eco. distress X closeness to M5S",
    "metropolitan_city" = "Metropolitan city"
  )


kable(cbind(names(var_dict_labels), var_dict_labels), row.names = FALSE,
      "latex", booktabs = T)  %>%
  kable_styling(font_size = 6)
```


# Introduction

```{r}
load("replication_file_01_ita_trust_unemployment.RData")
```

```{r ita_trust_unemployment, fig.width=9, fig.height=4}
grid.arrange(
  
  ggplot(ita_trust, aes(x=date, y=dont_trust_pp_perc)) + 
    scale_y_continuous(labels = percent) +
    geom_line() + 
    geom_point() + 
    geom_text(data = data.frame(x = as.Date("2018-07-15"), y=0.834, 
                                text = "Don't trust\npolitical parties"), 
              aes(x, y, label = text)) +
    scale_x_date(limits = c(min(ita_trust$date), as.Date("2018-12-15")),
                 breaks = as.Date(c("2000-01-01", "2005-01-01", 
                                    "2010-01-01", "2015-01-01")), 
                 date_labels = "%Y") +
    theme_bw() + 
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          panel.border = element_blank()) +
    labs(x=NULL, y=NULL),
  
  ggplot(ita_unemployment) + 
    geom_line(aes(x=TIME, y=Value,group=Territorio)) + 
    geom_point(aes(x=TIME, y=Value,group=Territorio)) + 
    scale_y_continuous(labels = percent) +
    geom_text(data = data.frame(x = as.Date("2018-07-15"),
                                y = c(0.11210570, 0.19350577),
                                text = c("Unemployment\n(Italy)",
                                         "Unemployment\n(Southern Italy)")),
              aes(x,y,label=text)) +
    scale_x_date(limits = c(min(ita_trust$date), as.Date("2018-12-15")),
                 breaks = as.Date(c("2000-01-01", "2005-01-01", 
                                    "2010-01-01", "2015-01-01")), 
                 date_labels = "%Y") +
    theme_bw() + 
    theme(legend.position = 'bottom',
          panel.border = element_blank()) + 
    labs(x=NULL, y=NULL),
  
  nrow = 2)
```

# Data

## Spatial analysis

Note: The following script is provided only to document the spatial analysis: **it is not replicable** given the data provided in this replication package. The output saved in the `replication_file_03_spatial_analysis.RData`. 

```{r, echo = TRUE}
read_chunk('replication_code_01_spatial_analysis.R')
```

```{r spatial-analysis, eval = FALSE}

```

```{r}
load('replication_file_03_spatial_analysis.RData')
```

## Meetup.com

```{r}
load('replication_file_02_meetup_event_venue.RData')

meetup_event$date <- 
  as.Date(meetup_event$time, "%Y-%m-%d %H:%M:%S")

# Unique number of groups
length(unique(meetup_group$group_id))
# Unique number of public groups 
length(unique(meetup_group$group_id[meetup_group$visibility == 'public']))
# Unique number of members
length(unique(meetup_group_member$member_id))
# Unique number of events
nrow(meetup_event)
# Unique number of venues
nrow(meetup_venue)
# Date of first event
min(meetup_event$date)
# Date of last event
max(meetup_event$date)
```

```{r}
kable(meetup_group %>% 
        dplyr::group_by(collected_on) %>%
        dplyr::summarize(n = n(),
                         `is public (%)` = sum(visibility == 'public') / n() * 100,
                         `n members (median)` = median(members),
                         `+50 members (%)` = sum(members >= 50) / n() * 100),
      digit = 2)

# 2014 groups missing in 2018
sum(!meetup_group$group_id[meetup_group$collected_on == '2014-08']
    %in% meetup_group$group_id[meetup_group$collected_on == '2018-07']) / 
  sum(meetup_group$collected_on == '2014-08') * 100
```

```{r meetup_group_visibility, fig.width = 13, fig.height = 8}
meetup_group_pnt.sp <- 
  SpatialPointsDataFrame(meetup_group[,c("lon","lat")], 
                         data = meetup_group[,c("lon","lat", "visibility", "collected_on",
                                                "members", "timestamp")],
                         proj4string = CRS(proj4string(comuni_sp)))
res <-
  over(meetup_group_pnt.sp, comuni_sp)

res <- 
  cbind(res, meetup_group[,c("visibility","collected_on", "timestamp","members")])

res$visibility <- as.factor(res$visibility)
levels(res$visibility)[levels(res$visibility)=="members"] <- "limited"
levels(res$visibility)[levels(res$visibility)=="public_limited"] <- "limited"

res$is_public <- res$visibility == "public"

res$north_south <- ifelse(as.numeric(as.character(res$COD_REG))<13, "North", "South")
res$north <- as.numeric(as.character(res$COD_REG))<13

res$dist_10k_ppl <- res$dist_10k_ppl + 1

res$days_from_origin <- as.numeric(as.Date(res$timestamp))

require(grid)
require(scales)

p1 <- 
  ggplot(res, aes(x = visibility, y = log(pop_density))) +
  geom_boxplot() +
  facet_grid(.~collected_on) + theme_bw() +
  labs(y = "Population density (log.)", x = "Meetup group visibility") +
  theme(axis.title.x = element_blank())
p2 <- 
  ggplot(res, aes(x = visibility, y = income2016_procapite)) +
  geom_boxplot() +
  facet_grid(.~collected_on) + theme_bw() +
  labs(y = "Income (procapite)", x = "Meetup group visibility")
p3 <- 
  res %>% 
  dplyr::filter(!is.na(north_south)) %>%
  dplyr::group_by(north_south, collected_on) %>% 
  dplyr::summarise(public_perc= sum(visibility == "public")/ n()) %>% 
  ggplot(aes(x = north_south, y = public_perc)) +
  geom_bar(stat='identity') +
  facet_grid(.~collected_on) + theme_bw() +
  labs(y = "Public groups (perc.)", x = NULL) +
  scale_y_continuous(labels = percent)
p4 <- 
  ggplot(res, aes(x = visibility, y = log(members))) +
  geom_boxplot() +
  facet_grid(.~collected_on) + theme_bw() +
  labs(y = "Members", x = "Meetup group visibility")

grid.newpage()
grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2), ggplotGrob(p3), ggplotGrob(p4), 
                size = "last"))
```

```{r}
dat <- res[!is.na(res$north),]
summary(dat$unemployment)
dat$dist_10k_ppl_log <- log(dat$dist_10k_ppl)
dat$members_log <- log(dat$members)

meetup_visibility_mod <-
  glm(is_public~members_log+
        dist_10k_ppl_log +
        unemployment+north+
        income2016_procapite+
        days_from_origin, 
      data = dat, 
      family = 'binomial', na.action="na.exclude")

newdata_north <- 
  data.frame(members_log = median(dat$members_log),
             dist_10k_ppl_log = median(dat$dist_10k_ppl_log), 
             unemployment = median(dat$unemployment),
             north = c(TRUE, FALSE),
             income2016_procapite = median(dat$income2016_procapite),
             days_from_origin = min(dat$days_from_origin))

newdata_unempl <- 
  data.frame(members_log = median(dat$members_log),
             dist_10k_ppl_log = median(dat$dist_10k_ppl_log), 
             unemployment = c(.05, .05, .2, .2),
             north = c(TRUE, FALSE, TRUE, FALSE),
             income2016_procapite = median(dat$income2016_procapite),
             days_from_origin = median(dat$days_from_origin))

newdata_dist_ppl <- 
  data.frame(members_log = median(dat$members_log),
             dist_10k_ppl_log = c(0.6931, 0.6931, 8.4502, 8.4502), 
             unemployment = median(dat$unemployment),
             north = c(TRUE, FALSE, TRUE, FALSE),
             income2016_procapite = median(dat$income2016_procapite),
             days_from_origin = median(dat$days_from_origin))

newdata_early <- 
  data.frame(members_log = median(dat$members_log),
             dist_10k_ppl_log = median(dat$dist_10k_ppl_log), 
             unemployment = c(.05, .05, .1, .1, .20, .20),
             north = c(TRUE, FALSE, TRUE, FALSE, TRUE, FALSE),
             income2016_procapite = median(dat$income2016_procapite),
             days_from_origin = min(dat$days_from_origin))

newdata_late <- 
  data.frame(members_log = median(dat$members_log),
             dist_10k_ppl_log = median(dat$dist_10k_ppl_log), 
             unemployment = c(.05, .05, .1, .1, .20, .20),
             north = c(TRUE, FALSE, TRUE, FALSE, TRUE, FALSE),
             income2016_procapite = median(dat$income2016_procapite),
             days_from_origin = max(dat$days_from_origin))

predict(meetup_visibility_mod, newdata_north, type="response")
predict(meetup_visibility_mod, newdata_unempl, type="response")
predict(meetup_visibility_mod, newdata_dist_ppl, type="response")
predict(meetup_visibility_mod, newdata_early, type="response")
predict(meetup_visibility_mod, newdata_late, type="response")

```

```{r, results = 'asis'}
do.call(
  stargazer, 
  c(list(meetup_visibility_mod,
          type = 'latex', out = "table/meetup_visibility.tex",
          dep.var.labels = "Meetup group is public",
          covariate.labels = c("Members (log.)", 
                               'Distance to first +10,000 pple. city (log.)',
                               'Unemployment rate (perc.)',
                               'North (binary)', 'Income (procapite)', 'Time (days)')),
    sg_param_appendix))
```


```{r, echo = TRUE}
read_chunk('replication_code_02_geocode_venue_italy.R')
```

```{r geocode-venue-italy, eval = FALSE}

```

```{r}
load('replication_file_04_geocode_venue_italy.RData')

kable(table(meetup_venue_it$regione))

```

```{r}

# Members 

meetup_member_by_data_collection <- data.frame()

for (d in unique(meetup_group_member$collected_on)) {
  this_meetup_member <- meetup_member[meetup_member$collected_on == d,]
  meetup_member_by_data_collection <- 
    rbind(meetup_member_by_data_collection, 
          data.frame(timestamp = d,
                     n = length(unique(this_meetup_member$member_id)),
                     missing = ifelse(exists('previous_meetup_member'),
                                      sum(!unique(this_meetup_member$member_id) %in%
                                            unique(previous_meetup_member$member_id)),
                                      NA)))
  previous_meetup_member <- this_meetup_member
}


meetup_group_member$joined <- as.Date(meetup_group_member$joined)

firstJoined <- function(x) {
  if (all(is.na(x))) {
    return(NA)
  } else {
    return(min(x))
  }
}

meetup_member_first_joined <-
  meetup_group_member %>%
  dplyr::group_by(member_id) %>%
  dplyr::summarise(first_joined = firstJoined(joined))

meetup_member_first_joined$first_joined[
  is.na(meetup_member_first_joined$first_joined)
  ] <-
  meetup_member$joined[
    match(meetup_member_first_joined$member_id[
      is.na(meetup_member_first_joined$first_joined)
      ] , meetup_member$member_id)]

meetup_member_by_month  <- 
  meetup_member_first_joined %>%
  dplyr::group_by(month = as.Date(cut(first_joined,  breaks = 'month'))) %>%
  dplyr::summarise(n = n())

meetup_member_by_month$n_cumsum <- cumsum(meetup_member_by_month$n)

meetup_member_by_data_collection$remaining <- 
  c(meetup_member_by_data_collection$n[1:4] - 
      meetup_member_by_data_collection$missing[2:5], NA)

meetup_member_by_month$type <- 'by first joining date available (cumulative)'
meetup_member_by_month <- 
  rbind(meetup_member_by_month, 
        data.frame(month = as.Date(paste0(meetup_member_by_data_collection$timestamp, "-15")),
                   n = NA,
                   n_cumsum = meetup_member_by_data_collection$n,
                   type = 'registered to existing groups\nat the time of data collection'))

meetup_event$time <- as.Date(meetup_event$time)
meetup_event_by_month <- 
  subset(meetup_event, yes_rsvp_count>0) %>%
  dplyr::group_by(month = as.Date(cut(time,  breaks = 'month'))) %>%
  dplyr::summarise(n = n())

plot_date_limits <- c(min(meetup_event_by_month$month)-90, 
                      max(meetup_event_by_month$month)-90)

# Events

meetup_event_by_month$n_ma12 <- 
  rollmean(meetup_event_by_month$n, k=12, fill=NA)

# Events/comune

meetup_venue_by_comune <- 
  merge(meetup_event[meetup_event$yes_rsvp_count>0,c('venue','time')],
        meetup_venue_it[,c('venue_id','istat_cod','regione')], 
        by.x = 'venue', by.y = 'venue_id')

meetup_venue_by_comune <- 
  meetup_venue_by_comune %>%
  dplyr::group_by(month = cut(time, breaks = 'month'),
                  istat_cod,
                  regione) %>%
  dplyr::summarise(n = n())

meetup_venue_by_comune$pop2011 <- NA
meetup_venue_by_comune$pop2011 <- 
  comuni_sp$pop2011[match(meetup_venue_by_comune$istat_cod, 
                          comuni_sp$PRO_COM_T)]

meetup_venue_by_month <-
  meetup_venue_by_comune %>%
  dplyr::group_by(month) %>%
  dplyr::summarize(perc_pop = 
                     (sum(pop2011,na.rm=T)/sum(comuni_sp$pop2011)))
meetup_venue_by_month$month <- 
  as.Date(as.character(meetup_venue_by_month$month))

meetup_venue_by_month$perc_pop_ma12 <- 
  rollmean(meetup_venue_by_month$perc_pop, k=12,fill=NA)



p1 <- 
  ggplot(meetup_event_by_month, 
         aes(x=month,y=n)) +
  geom_line() +
  geom_line(aes(y=n_ma12), linetype = 2, color = 'gray') + 
  theme_bw() + labs(y=NULL, x='Meetup monthly events') +
  scale_x_date(limits = plot_date_limits, date_breaks = 'year') + 
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank())

p2 <-
  ggplot(meetup_member_by_month, aes(x=month, y=n_cumsum)) +
  geom_line(aes(linetype = type)) + scale_y_continuous(labels = comma) +
  geom_point(data = meetup_member_by_data_collection, 
             aes(x=as.Date(paste0(timestamp, "-15")), y=n), stat ='identity') +
  theme_bw() + labs(y=NULL, x='Meetup users') +
  scale_x_date(limits = plot_date_limits, date_breaks = 'year', date_labels = "%Y") +
  theme(legend.title=element_blank(),
        legend.position = c(0.2, 0.8),
        legend.background = element_blank())

p3 <- 
  ggplot(subset(meetup_venue_by_month, month < plot_date_limits[2] - 30), 
         aes(x=month, y=perc_pop)) +
  geom_line() +
  geom_line(aes(y=perc_pop_ma12), linetype = 2, color = 'gray') + 
  theme_bw() + 
  labs(y=NULL, 
       x='Population in comune with one monthly event or more (over total population)') +
  scale_x_date(limits = plot_date_limits, date_breaks = 'year') +  
  scale_y_continuous(labels = percent) + 
  theme(axis.ticks.x = element_blank(),
        legend.title=element_blank(),
        legend.position = c(0.1, 0.7),
        axis.text.x = element_blank(),
        legend.background = element_blank())


```

```{r}
read_chunk('replication_code_05_meetup_frequency.R')
```

```{r meetup-frequency, eval = FALSE}

```

```{r meetup_7days_comune_8w_win, width = 10, height = 5.5}

load("06_meetup_7days_comune_8w_win.RData")
meetup_7days_coeff$p.value <- 
  as.numeric(cut(meetup_7days_coeff$meetup_7days_pvalue, 
                 breaks = c(0,.01,.02,.03,.04,.05,.1,1), labels = 1:7))

ggplot(meetup_7days_coeff, aes(x=date)) +
  geom_line(aes(y=4.8, color=p.value), size=8) +
  scale_color_gradientn(colours = viridis(12, direction=-1)[c(1:7,12)],
                        labels = c('0', '.01', '.02', 
                                   '.03', '.04', '.05','.1')) +
  geom_line(aes(y=meetup_7days_estimate*100)) +
  geom_ribbon(aes(ymin=lower*100,ymax=upper*100), alpha=.3) +
  geom_vline(xintercept = as.Date(c("2013-02-24","2018-03-04"))) + 
  geom_vline(xintercept = as.Date(
    c("2012-05-07", "2012-10-28", "2014-05-25")), 
    linetype = 2, alpha = .3) + 
  labs(y=NULL, x="Effect on M5S vote in the following election") +
  facet_grid(north_south ~ .) + 
  theme_bw() + theme(legend.position = 'bottom')
```

### Figure 2

```{r meetup_timeseries, fig.width=9, fig.height=5.5}
plot_grid(p3, p1, p2, align = "v", nrow = 3, rel_heights = c(1/3, 1/3, 1/3))
```

```{r m5s_2018_2013_compar_meetup, fig.width=20, fig.height=15}
load('replication_file_05_regione2018_sp_sf.RData')

hex_grid_25km_meetup2014$Freq[is.na(hex_grid_25km_meetup2014$Freq)] <- 0
hex_grid_25km_meetup2018$Freq[is.na(hex_grid_25km_meetup2018$Freq)] <- 0

hex_grid_25km_meetup2014$brks2014 <- 
  cut(hex_grid_25km_meetup2014$Freq,
      breaks = c(-1, 0, 1, 5, 10, 300))
hex_grid_25km_meetup2018$brks2018 <- 
  cut(hex_grid_25km_meetup2018$Freq,
      breaks = c(-1, 0, 1, 5, 10, 300))

hex_grid_25km_meetup2014_sf <- st_as_sf(hex_grid_25km_meetup2014)
hex_grid_25km_meetup2018_sf <- st_as_sf(hex_grid_25km_meetup2018)

grid.arrange(
  ggplot() +
    geom_sf(data=hex_grid_25km_meetup2014_sf, aes(fill=brks2014), color = NA) +
    scale_fill_brewer(palette='BuGn') +
    geom_sf(data = regione2018_sf, colour = 'black', fill = NA, alpha = .7) +
    geom_point(data=as.data.frame(meetup2014_event), aes(lon,lat), size = .2, alpha = .8) +
    coord_sf() +
    blank() +
    guides(fill=guide_legend(title="Meetup event count")) +
    labs(title = '2013') +
    theme(legend.position="bottom") +
    theme(legend.key.width=unit(2, "cm")),
  ggplot() +
    geom_sf(data=hex_grid_25km_meetup2018_sf, aes(fill=brks2018), color = NA) +
    scale_fill_brewer(palette='BuGn') +
    geom_sf(data = regione2018_sf, colour = 'black', fill = NA, alpha = .7) +
    geom_point(data=as.data.frame(meetup2018_event), aes(lon,lat), size = .2, alpha = .8) +
    coord_sf() +
    blank() +
    guides(fill=guide_legend(title="Meetup event count")) +
    labs(title = '2018') +
    theme(legend.position="bottom") +
    theme(legend.key.width=unit(2, "cm")),
  ncol = 2
)
```

## Election data

```{r}
load("replication_file_06_regione2018_south_sp.RData")
istat_rip_geo_2018_sf <- st_as_sf(istat_rip_geo_2018)
st_crs(istat_rip_geo_2018_sf) <- 32632

comuni_sf_keep005 <- st_as_sf(comuni_sp_keep005)

comuni_sf_keep005 <- merge(comuni_sf_keep005, comuni_sp, by = 'PRO_COM_T')

comuni_sp$M5S_diff <- comuni_sp$M5S_diff / 100

## Compare 2013 and 2018 
### Difference from national mean

cut_breaks <- 
  c(-0.5, -0.3, -0.1, -0.05, 0.05, 0.1, 0.3, 0.5)
cut_labels <- 
  c("-50%,-30%", "-30%,-10%", "-10%,-5%","-5%,+5%", 
    "+5%,+10%","+10%,+30%","+30%,+50%")
nat_vote_2013 <- 
  0.2566
nat_vote_2018 <- 
  0.3267

plot2013 <-
  ggplot() + 
  geom_sf(data=comuni_sf_keep005,
          aes(fill =
                cut(
                  M5S_2013 -
                    nat_vote_2013,
                  breaks = cut_breaks,
                  labels = cut_labels)),
          col = NA) +
  geom_sf(data = regione2018_sf, 
          colour = 'gray', fill = NA,alpha=.5) +
  geom_sf(data = istat_rip_geo_2018_sf[istat_rip_geo_2018_sf$FID>2,],
          colour = 'black', fill = NA, size=1) +
  scale_fill_brewer(palette = 'Spectral', drop = FALSE, direction=-1) +
  coord_sf() +
  blank() +
  guides(fill=FALSE) +
  theme(plot.caption = element_text(hjust=0.5, size=21)) +
  labs(caption = paste0('2013 (national vote: ', 
                        round(nat_vote_2013*100,1),'%)'))

plot2018 <-
  ggplot() + 
  geom_sf(data=comuni_sf_keep005, 
          aes(fill = 
                cut(
                  M5S_2018 - 
                    nat_vote_2018,
                  breaks = cut_breaks,
                  labels = cut_labels)), 
          col = NA) +
  geom_sf(data = regione2018_sf, colour = 'gray', fill = NA,alpha=.5) +
  scale_fill_brewer(palette = 'Spectral', drop = FALSE, direction=-1) +
  coord_sf() +
  blank() +
  guides(fill=guide_legend(nrow=1,byrow=TRUE,title="Difference from national vote",
                           label.position = "bottom",
                           label.hjust = 5)) +
  theme(legend.position="bottom",legend.direction="horizontal",
        legend.title=element_blank(),
        legend.text = element_text(size = 18, margin = margin(t = 5)),
        plot.title = element_text(size = 30),
        legend.spacing.x = unit(2, 'cm'),
        legend.key = element_rect(colour = 'black', fill = NA, size = 0.5)) +
  labs(title = paste0('2018 (national vote: ', round(nat_vote_2018*100,2),'%)'))

extractLegend <- function(a.gplot) {
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}
plot_legend <- extractLegend(plot2018)

plot2018 <-
  ggplot() + 
  geom_sf(data=comuni_sf_keep005,
          aes(fill =
                cut(
                  M5S_2018 -
                    nat_vote_2018,
                  breaks = cut_breaks,
                  labels = cut_labels)),
          col = NA) +
  geom_sf(data = regione2018_sf, colour = 'gray', fill = NA, alpha=.5) +
  geom_sf(data = istat_rip_geo_2018_sf[istat_rip_geo_2018_sf$FID>2,],
          colour = 'black', fill = NA, size=1) +
  scale_fill_brewer(palette = 'Spectral', drop = FALSE, direction=-1) +
  coord_sf() +
  blank() +
  guides(fill=FALSE) +
  theme(plot.caption = element_text(hjust=0.5, size=21)) +
  labs(caption = paste0('2018 (national vote: ', round(nat_vote_2018*100,2),'%)'))

```

### Figure 1 (Supplemental online material)

```{r m5s_2018_2013_diff_from_nat_mean,  fig.width=12, fig.height=8}
grid.arrange(
  plot2013,
  plot2018,
  plot_legend,
  heights = c(8, 1),
  layout_matrix = rbind(c(1,2),c(3))
)
```

```{r m5s_2018_diff_from_2013, fig.width=9, fig.height=12}
ggplot() + 
  geom_sf(data=comuni_sf_keep005,
          aes(fill = cut(M5S_2018-M5S_2013,
                         breaks = c(-.5, -.1, -.05, .05, 1, .5),
                         labels = c("-50%,-10%", "-10%,-5%", "-5%,+5%", "+5%,+10%", "+10%,+50%"),
                         include.lowest = T)),
          col = NA) +
  scale_fill_brewer(type="seq", palette = "PuOr", drop = FALSE) +
  geom_sf(data = regione2018_sf, fill = NA, alpha=.7) +
  geom_sf(data = istat_rip_geo_2018_sf[istat_rip_geo_2018_sf$FID>2,],
          colour = 'black', fill = NA, size=1) +
  coord_sf() +
  blank() +
  guides(fill = guide_legend(label.position = "bottom")) +
  theme(legend.position="bottom",
        legend.direction="horizontal",
        legend.title=element_blank(),
        legend.text = element_text(size = 18, margin = margin(t = 5)),
        legend.spacing.x = unit(.5, 'cm'),
        legend.key = element_rect(colour = 'black', fill = NA, size = 0.5))
```

## Opinion survey data

```{r, echo = TRUE}
read_chunk('replication_code_03_prov_demographics.R')
```

```{r prov-demographics, eval = FALSE}

```

```{r, echo = TRUE}
read_chunk('replication_code_04_reg_demographics.R')
```

```{r reg-demographics, eval = FALSE}

```

```{r income_unemployment_nuts2, fig.width = 5, fig.height = 5}
dat1 <- 
  eurostat::get_eurostat("tgs00026", 
                         time_format = "raw", stringsAsFactors = FALSE)
dat2 <- 
  eurostat::get_eurostat("lfst_r_lfu3rt", 
                         time_format = "raw", stringsAsFactors = FALSE) 

datit <- 
  merge(dat1 %>% dplyr::filter(nchar(geo) == 4, grepl("^IT", geo)), 
        dat2 %>% dplyr::filter(age == 'Y15-74', 
                               sex == 'T', nchar(geo) == 4, grepl("^IT", geo)),
        by = c('geo','time'))
datit$northsouth <- 
  ifelse(grepl("^ITF|^ITG", datit$geo), "South", "North")

datit$south <- 0
datit$south[datit$northsouth == 'South'] <- 1

dat1 <- 
  dat1 %>% dplyr::filter(nchar(geo) == 4, time == 2015)

dat2 <- 
  dat2 %>% dplyr::filter(age == 'Y15-74', sex == 'T', 
                         time == 2015, nchar(geo) == 4)

table(dat1$time)

merge_dat <- merge(dat1, dat2, by = 'geo')

merge_dat$it <- grepl("^IT", merge_dat$geo)
merge_dat$`Italian region` <- NA
merge_dat$`Italian region`[grepl("^ITF|^ITG", merge_dat$geo)] <- 'South'
merge_dat$`Italian region`[
  grepl("^ITC|^ITD|^ITE|^ITI|^ITH", merge_dat$geo)] <- 
  'North'

merge_dat$cnt <- str_extract(merge_dat$geo, "^[A-Z]{2}")

p <- 
  ggplot() +
  geom_point(data=merge_dat[merge_dat$it,], 
             aes(x=values.x,y=values.y/100, 
                 shape = `Italian region`), size = 2) +
  geom_point(data=merge_dat[!merge_dat$it,], 
             aes(x=values.x,y=values.y/100,fill=it), 
             alpha = 0.6, size = .6, shape = 3) + 
  scale_colour_manual(values = c('red','blue')) +
  scale_fill_manual(labels = NULL, values = 'black') +
  theme_bw() +
  scale_x_continuous(limits = c(5000, 25000), label = comma) +
  scale_y_continuous(label = percent) +
  labs(x = 'Disposable household income', y = 'Unemployment', fill = 'non-Italian region') + 
  theme(legend.position = 'bottom')

ggExtra::ggMarginal(p, type = "boxplot")
```

### Summary statistics of Istat surveys

```{r}
load("secondary_data/istat_aspetti_survey/AVQ_2013_IT_R/MICRODATI/DF_AVQ_A2013.RData")
load("secondary_data/istat_aspetti_survey/AVQ_2014_IT_R/MICRODATI/DF_AVQ_A2014.RData")
load("secondary_data/istat_aspetti_survey/AVQ_2015_IT_R/MICRODATI/DF_AVQ_A2015.RData")
load("secondary_data/istat_aspetti_survey/AVQ_2016_IT_R/MICRODATI/DF_AVQ_A2016.RData")

DF_AVQ_A2013$nat_pol_trust <- 
  rescale(DF_AVQ_A2013$puntifi4 + DF_AVQ_A2013$puntifi1, 
          to = c(0,10),
          from = c(0,20))
DF_AVQ_A2013$loc_pol_trust <- 
  rescale(DF_AVQ_A2013$puntifi8 + DF_AVQ_A2013$puntifi9 +
            DF_AVQ_A2013$puntifi10,
          to = c(0,10),
          from = c(0,30))
DF_AVQ_A2013$pol_trust <- 
  with(DF_AVQ_A2013, 
       rescale(puntifi8 + puntifi9 + puntifi10 + 
                 puntifi4 + puntifi1,
               to = c(0,10),
               from = c(0,50)))

DF_AVQ_A2014$nat_pol_trust <- 
  rescale(DF_AVQ_A2014$PUNTIFI4 + DF_AVQ_A2014$PUNTIFI1, 
          to = c(0,10),
          from = c(0,20))
DF_AVQ_A2014$loc_pol_trust <- 
  rescale(DF_AVQ_A2014$PUNTIFI8 + DF_AVQ_A2014$PUNTIFI9 + 
            DF_AVQ_A2014$PUNTIFI10,
          to = c(0,10),
          from = c(0,30))
DF_AVQ_A2014$pol_trust <- 
  with(DF_AVQ_A2014, 
       rescale(PUNTIFI8 + PUNTIFI9 + PUNTIFI10 + 
                 PUNTIFI4 + PUNTIFI1,
               to = c(0,10),
               from = c(0,50)))

DF_AVQ_A2015$nat_pol_trust <- 
  rescale(DF_AVQ_A2015$PUNTIFI4 + DF_AVQ_A2015$PUNTIFI1, 
          to = c(0,10),
          from = c(0,20))
DF_AVQ_A2015$loc_pol_trust <- 
  rescale(DF_AVQ_A2015$PUNTIFI8 + DF_AVQ_A2015$PUNTIFI9 + 
            DF_AVQ_A2015$PUNTIFI10,
          to = c(0,10),
          from = c(0,30))
DF_AVQ_A2015$pol_trust <- 
  with(DF_AVQ_A2015, 
       rescale(PUNTIFI8 + PUNTIFI9 + PUNTIFI10 + 
                 PUNTIFI4 + PUNTIFI1,
               to = c(0,10),
               from = c(0,50)))

DF_AVQ_A2016$nat_pol_trust <- 
  rescale(DF_AVQ_A2016$PUNTIFI4 + DF_AVQ_A2016$PUNTIFI1, 
          to = c(0,10),
          from = c(0,20))
DF_AVQ_A2016$loc_pol_trust <- 
  rescale(DF_AVQ_A2016$PUNTIFI8 + DF_AVQ_A2016$PUNTIFI9 + 
            DF_AVQ_A2016$PUNTIFI10,
          to = c(0,10),
          from = c(0,30))
DF_AVQ_A2016$pol_trust <- 
  with(DF_AVQ_A2016, 
       rescale(PUNTIFI8 + PUNTIFI9 + PUNTIFI10 + 
                 PUNTIFI4 + PUNTIFI1,
               to = c(0,10),
               from = c(0,50)))

## 2013
weighted.mean(DF_AVQ_A2013$nat_pol_trust, 
              weight = DF_AVQ_A2013$COEFIN, na.rm = T)
weighted.mean(DF_AVQ_A2013$loc_pol_trust, 
              weight = DF_AVQ_A2013$COEFIN, na.rm = T)
weighted.mean(DF_AVQ_A2013$pol_trust, 
              weight = DF_AVQ_A2013$COEFIN, na.rm = T)

## 2014
weighted.mean(DF_AVQ_A2014$nat_pol_trust, 
              weight = DF_AVQ_A2014$COEFIN, na.rm = T)
weighted.mean(DF_AVQ_A2014$loc_pol_trust, 
              weight = DF_AVQ_A2014$COEFIN, na.rm = T)
weighted.mean(DF_AVQ_A2014$pol_trust, 
              weight = DF_AVQ_A2014$COEFIN, na.rm = T)

## 2015
weighted.mean(DF_AVQ_A2015$nat_pol_trust, 
              weight = DF_AVQ_A2015$COEFIN, na.rm = T)
weighted.mean(DF_AVQ_A2015$loc_pol_trust, 
              weight = DF_AVQ_A2015$COEFIN, na.rm = T)
weighted.mean(DF_AVQ_A2015$pol_trust, 
              weight = DF_AVQ_A2015$COEFIN, na.rm = T)

## 2016
weighted.mean(DF_AVQ_A2016$nat_pol_trust, 
              weight = DF_AVQ_A2016$COEFIN, na.rm = T)
weighted.mean(DF_AVQ_A2016$loc_pol_trust, 
              weight = DF_AVQ_A2016$COEFIN, na.rm = T)
weighted.mean(DF_AVQ_A2016$pol_trust, 
              weight = DF_AVQ_A2016$COEFIN, na.rm = T)

```

```{r}
load("secondary_data/istat_aspetti_survey/AVQ_2013_IT_R/MICRODATI/DF_AVQ_A2013.RData")

DF_AVQ_A2013$reg[DF_AVQ_A2013$reg %in% c(41,42)] <- 300

DF_AVQ_A2013$coemicro <- 
  as.numeric(as.character(DF_AVQ_A2013$coemicro))
# weighted.mean(as.numeric(dat$ppapo==2), weight=coemicro, na.rm=T)

# Correlation local / national
with(DF_AVQ_A2013, cor.test(puntifi4 + puntifi1, puntifi8 + puntifi9 + puntifi10))
DF_AVQ_A2013$pol_trust_diff <- 
  with(DF_AVQ_A2013, (puntifi1 + puntifi4) - (puntifi8 + puntifi9 + puntifi10))

dat <- DF_AVQ_A2013

# Political trust
dat$loc_pol_trust <- rescale(with(dat, (puntifi8 + puntifi9 + puntifi10)), 
                             to = c(0,10),
                             from = c(0,30))
dat$nat_pol_trust <- rescale(with(dat, (puntifi1 + puntifi4),
                                  to = c(0,10),
                                  from = c(0,20)))

asIndexComp <- function(x) {
  x <- as.numeric(x)
  x[is.na(x)] <- 0
  x
}
reverse <- function(x) {
  if (is.na(x)) {
    res <- NA 
  } else if (x == 1) {
    res <- 3
  } else if (x == 2) {
    res <- 2
  } else if (x == 3) {
    res <- 1
  } else if (x == 4) {
    res <- 0
  } else {
    res <- NA
  }
  return(res)
}

# Social capital: association membership

dat$member <-
  as.logical(as.numeric(with(dat,
                             asIndexComp(volon==4) + 
                               asIndexComp(psind == 4) + asIndexComp(pgrvo == 6) + 
                               asIndexComp(paeco == 2) + asIndexComp(pcult == 4) +
                               asIndexComp(paspro == 6) > 0)))

# Social capital: Trust  
dat$interpers_trust <-
  rescale(
    with(dat,
         asIndexComp(fiducia == 1) + 
           sapply(fidu1, reverse)*0.1 + 
           sapply(fidu2, reverse)*0.2 + 
           sapply(fidu3, reverse)*0.4),
    to = c(0, 10),
    from = c(0, 3.1))


dat$edu <- 
  as.numeric(factor(as.factor(dat$istr), 
                    levels = c('11', '10', '9', '7', '2', '1')))
dat$north <- dat$rip < 4
dat$unemployed <- dat$cond %in% 2:3
dat$type_comune <- factor(as.character(dat$dom))
levels(dat$type_comune) <- c('Metropolitan area', 'Periphery of metropolitan area', 
                             'Comune up to 2000', 'Comune between 2001 and 10,000',
                             'Comune between 10.001 and 50,000', 'Comune over 50,000')
dat$male <- dat$sesso == 1
dat$age <- dat$eta

# Political interest
dat$pol_interest <-
  rescale(with(dat, polit), to = c(0,10), from = c(6,1))

# Political participation
dat$pol_part <- 
  rescale(
    with(dat, 
       as.numeric(comiz==2) + 
         as.numeric(cortei==4) + 
         as.numeric(dibpo==6) +
          as.numeric(ppapo==2) +
          as.numeric(volpa==8)),
    to = c(0,10))


aspetti2013_ind_lev_lmmod1 <- 
  lm(data = dat, pol_part ~ nat_pol_trust + loc_pol_trust + 
       edu + north + member + 
       interpers_trust + 
       type_comune + unemployed + age + male)

aspetti2013_ind_lev_lmmod3 <- 
  lm(data = dat, loc_pol_trust ~ nat_pol_trust + edu + member + north + 
       type_comune + unemployed + age + male)

```

### Table 7 (Supplemental online material)

```{r, results = 'asis'}

do.call(stargazer,
        c(list(aspetti2013_ind_lev_lmmod1, aspetti2013_ind_lev_lmmod3,
               type = 'latex', out = "table/aspetti2013_ind_lev_lmmod_full.tex",
               dep.var.labels = var_dict_labels[c("pol_part", "loc_pol_trust")],
               column.labels = c("(B)", "(B)"),
               covariate.labels = var_dict_labels[
                 c('nat_pol_trust', 'loc_pol_trust',
                   'edu', 'north',
                   'member', 'interpers_trust',
                   "type_comunePeriphery of metropolitan area",
                   "type_comuneComune up to 2000",
                   "type_comuneComune between 2001 and 10,000",
                   "type_comuneComune between 10.001 and 50,000",
                   "type_comuneComune over 50,000",
                   'unemployed', 'age', 'male')]),
          sg_param_appendix))
```

```{r}
load('replication_file_08_region_sp_list.RData')

region_2013_sf <- st_as_sf(region_sp_list[[1]])
region_2014_sf <- st_as_sf(region_sp_list[[2]])
region_2015_sf <- st_as_sf(region_sp_list[[3]])
region_2016_sf <- st_as_sf(region_sp_list[[4]])

region_2013_sf$nat_pol_trust <- 
  rescale(region_2013_sf$puntifi4_mean + region_2013_sf$puntifi1_mean, 
          to = c(0,10),
          from = c(0,20))
region_2013_sf$loc_pol_trust <- 
  rescale(region_2013_sf$puntifi8_mean + region_2013_sf$puntifi9_mean +
            region_2013_sf$puntifi10_mean,
          to = c(0,10),
          from = c(0,30))
region_2013_sf$pol_trust <- 
  with(region_2013_sf, 
       rescale(puntifi8_mean + puntifi9_mean + puntifi10_mean + 
                 puntifi4_mean + puntifi1_mean,
               to = c(0,10),
               from = c(0,50)))

region_2014_sf$nat_pol_trust <- 
  rescale(region_2014_sf$puntifi4_mean + region_2014_sf$puntifi1_mean, 
          to = c(0,10),
          from = c(0,20))
region_2014_sf$loc_pol_trust <- 
  rescale(region_2014_sf$puntifi8_mean + region_2014_sf$puntifi9_mean + 
            region_2014_sf$puntifi10_mean,
          to = c(0,10),
          from = c(0,30))
region_2014_sf$pol_trust <- 
  with(region_2014_sf, 
       rescale(puntifi8_mean + puntifi9_mean + puntifi10_mean + 
                 puntifi4_mean + puntifi1_mean,
               to = c(0,10),
               from = c(0,50)))

region_2015_sf$nat_pol_trust <- 
  rescale(region_2015_sf$puntifi4_mean + region_2015_sf$puntifi1_mean, 
          to = c(0,10),
          from = c(0,20))
region_2015_sf$loc_pol_trust <- 
  rescale(region_2015_sf$puntifi8_mean + region_2015_sf$puntifi9_mean + 
            region_2015_sf$puntifi10_mean,
          to = c(0,10),
          from = c(0,30))
region_2015_sf$pol_trust <- 
  with(region_2015_sf, 
       rescale(puntifi8_mean + puntifi9_mean + puntifi10_mean + 
                 puntifi4_mean + puntifi1_mean,
               to = c(0,10),
               from = c(0,50)))

region_2016_sf$nat_pol_trust <- 
  rescale(region_2016_sf$puntifi4_mean + region_2016_sf$puntifi1_mean, 
          to = c(0,10),
          from = c(0,20))
region_2016_sf$loc_pol_trust <- 
  rescale(region_2016_sf$puntifi8_mean + region_2016_sf$puntifi9_mean + 
            region_2016_sf$puntifi10_mean,
          to = c(0,10),
          from = c(0,30))
region_2016_sf$pol_trust <- 
  with(region_2016_sf, 
       rescale(puntifi8_mean + puntifi9_mean + puntifi10_mean + 
                 puntifi4_mean + puntifi1_mean,
               to = c(0,10),
               from = c(0,50)))

is <- c('4','1','8','9','10')
labels <- c('1' = 'Parliament', '4' = 'Parties', '8' = 'Regional gov.', 
            '9' = 'Provincial gov.', '10' = 'Municipal gov.')
```

```{r, fig.show = 'hide'}
plotRegIndex <- function(year) {
  this_sf <- get(sprintf('region_%s_sf', year))
  this_sf$n_meetup_events_per_10k <- this_sf$n_meetup_events / (this_sf$population/10000)
  this_mean_nat <- round(mean(this_sf$nat_pol_trust),2)
  this_mean_loc <- round(mean(this_sf$loc_pol_trust),2)
  this_mean_trust <- round(mean(this_sf$pol_trust),2)
  this_mean_meetups <- round(mean(this_sf$n_meetup_events_per_10k),2)
  grid.arrange(
    ggplot(this_sf) + 
      geom_sf(aes(fill = nat_pol_trust), colour = 'black', size = .04) +
      scale_fill_gradientn(colours = brewer.pal(7, "PuOr")) +
      theme_bw() +
      labs(caption = paste0("Trust in national inst.", 
                            " (avg.=",this_mean_nat,")"), 
           fill = NULL, title = year) +
      theme(axis.text = element_blank(),
            axis.ticks = element_blank(), 
            legend.position = c(0.5, 0.02),
            legend.direction="horizontal",
            legend.key.width=unit(.8,"cm"),
            legend.key.height=unit(.15,"cm")),
    ggplot(this_sf) + 
      geom_sf(aes(fill = loc_pol_trust), colour = 'black', size = .04) +
      scale_fill_gradientn(colours = brewer.pal(7, "PuOr")) +
      theme_bw() +
      labs(caption = paste0("Trust in local inst.", 
                            " (avg.=",this_mean_loc,")"), 
           fill = NULL, title = NULL) +
      theme(axis.text = element_blank(),
            axis.ticks = element_blank(), 
            legend.position = c(0.5, 0.02),
            legend.direction="horizontal",
            legend.key.width=unit(.8,"cm"),
            legend.key.height=unit(.15,"cm")),
    ggplot(this_sf) + 
      geom_sf(aes(fill = pol_trust), colour = 'black', size = .04) +
      scale_fill_gradientn(colours = brewer.pal(7, "PuOr")) +
      theme_bw() +
      labs(caption = paste0("Overal trust in inst.", 
                            " (avg.=",this_mean_trust,")"), 
           fill = NULL, title = NULL) +
      theme(axis.text = element_blank(),
            axis.ticks = element_blank(), 
            legend.position = c(0.5, 0.02),
            legend.direction="horizontal",
            legend.key.width=unit(.8,"cm"),
            legend.key.height=unit(.15,"cm")),
    ggplot(this_sf) + 
      geom_sf(aes(fill = n_meetup_events_per_10k), 
              colour = 'black', size = .04) +
      scale_fill_gradientn(colours = brewer.pal(7, "YlGnBu")) +
      theme_bw() +
      labs(caption = paste0("Events per 10k people", 
                            " (avg.=",this_mean_meetups,")"), 
           fill = NULL, title = NULL) +
      theme(axis.text = element_blank(),
            axis.ticks = element_blank(), 
            legend.position = c(0.5, 0.02),
            legend.direction="horizontal",
            legend.key.width=unit(.8,"cm"),
            legend.key.height=unit(.15,"cm")),
    nrow = 4
  )
}

plots <- lapply(2013:2016, plotRegIndex)
```

### Figure 11 (Supplemental online material)

```{r ita_trust_index_reg_ts_map, fig.width = 9, fig.height = 9}
grid.arrange(grobs = plots, ncol = 4)
```

```{r onlinesurvey_issues, fig.width =  8.5, fig.height = 7, eval = FALSE}
dat <- read_excel("secondary_data/bfna/Disruptive_Italy 2018.xlsx", 
                  sheet = "baza date completes")
dat <- dat[!is.na(dat$Q1),]

# Coding
dat$m5s_voter <- dat$Q42 == 1
dat$SouthIslands <- dat$Q3 %in% c("Islands","South")

dat$age <- as.numeric(as.Date('2018-03-04') - as.Date(dat$Q2)) / 365

dat$Q14[is.na(dat$Q14)] <- 7
dat$Q14[dat$Q14 == 997] <- 7

dat$Q24[dat$Q24 == 997] <- 7

dat$Q32_3[dat$Q32_3 == 999] <- 1

dat$Q33[dat$Q33] <- dat$Q33 < 3

dat$Q34[dat$Q34] <- dat$Q34 < 3

dat$Q39[dat$Q39 == 999] <- NA


dat$Q46 <- dat$Q46 < 3

dat$Q48 <- dat$Q48 < 3

dat$Q54_8[dat$Q54_8 == 999] <- NA

ivars <- c('Q1','age','Q5','Q6', 'Q14', 'Q24',
           paste0('Q27_bis_', 1:6),
           paste0('Q28_0', 1:8),
           paste0('Q29_', 1:9),
           paste0('Q30_0', 1:10),
           paste0('Q31_2_0', 1:9),
           paste0('Q32_0', 1:11),
           'Q32_1', 'Q32_3',
           'Q33', 'Q34', 'Q36', 
           paste0('Q38_0', 1:11),
           paste0('Q45_0', 1:10),
           'Q46', 'Q47', 'Q48',
           paste0('Q53_', 1:9),
           paste0('Q54_', 1:8))

for (var in ivars) {
  if (length(table(dat[[var]])) == 1) {
    dat[[var]] <- dat[[var]] %in% 1
  }
}

q_demo <- 1:6
q_mediause <- 11:29
q_civicpart <- 30:54

dat_m5s <- subset(dat, m5s_voter)

# Significance

require(ggplot2)
require(scales)

dat$SouthIslands <- as.factor(dat$SouthIslands)
levels(dat$SouthIslands) <- c("North", "South")

dat$m5s_voter <- as.character(dat$m5s_voter)
dat$m5s_voter[is.na(dat$m5s_voter)] <- "Didn't\nvote"
dat$m5s_voter <- as.factor(dat$m5s_voter)
levels(dat$m5s_voter)[2:3] <- c("Other\nparty", "M5S")
dat$m5s_voter <- factor(dat$m5s_voter, levels = rev(levels(dat$m5s_voter)))

## All issues
require(reshape2)
require(dplyr)
dat_issues <- dat[,grepl("^Q32_0\\d{1,2}$|m5s_voter|SouthIslands", names(dat))]
dat_issues <- melt(dat_issues, id.vars = c("m5s_voter","SouthIslands"))
dat_issues <-
  dat_issues %>% 
  dplyr::group_by(variable) %>%
  dplyr::mutate(nat_true = sum(value)/n()) %>%
  dplyr::group_by(m5s_voter, SouthIslands, variable) %>%
  dplyr::summarise(group_diff =  sum(value)/n() - nat_true[1],
                   group_sd   = sd(value),
                   se   = sd(value) / sqrt(n()))

levels(dat_issues$variable) <- 
  c('Corruption', 'Costs of living', 'Discrimination', 'Education', 
    'Healthcare', 'Lack of security', 'Lack of economic stability',
    'Immigration', 'Political instability', 'Terrorism', 'Unemployment')

ggplot(dat_issues, 
       aes(fill=SouthIslands, y=group_diff, x=m5s_voter)) +
  geom_bar(stat = 'identity', position = 'dodge') +
  geom_errorbar(aes(ymin=group_diff-se, ymax=group_diff+se),
                width=.2, position=position_dodge(.9)) +
  scale_y_continuous(labels = scales::percent) +
  facet_wrap(~variable) + theme_bw() +
  theme(legend.position = 'bottom',
        legend.title=element_blank()) +
  labs(x='In 2018 voted...', y=NULL)
```

## Models

```{r}
camera_sp$north <- 
  camera_sp$macro_region %in% 
  c('North-west', 'North-east', 'Centre')
comuni_sp$north <- 
  comuni_sp$macro_region %in% 
  c('North-west', 'North-east', 'Centre')

## 2018
camera_sp$M5S_2018_neigh_avg <- NA
for (i in 1:nrow(camera_sp)) {
  camera_sp$M5S_2018_neigh_avg[i] <- 
    mean(camera_sp$M5S_2018[camera_nb[['10']][[i]]])
}
camera_sp$M5S_2018_neigh_avg[
  is.na(camera_sp$M5S_2018_neigh_avg)] <- 0 # Islands

comuni_sp$M5S_2018_neigh_avg <- NA
for (i in 1:nrow(comuni_sp)) {
  comuni_sp$M5S_2018_neigh_avg[i] <- 
    mean(comuni_sp$M5S_2018[comuni_nb[['10']][[i]]])
}
comuni_sp$M5S_2018_neigh_avg[
  is.na(comuni_sp$M5S_2018_neigh_avg)] <- 0 # Islands

# 2013
comuni_sp$M5S_2013_neigh_avg <- NA
for (i in 1:nrow(comuni_sp)) {
  comuni_sp$M5S_2013_neigh_avg[i] <- 
    mean(comuni_sp$M5S_2013[comuni_nb[['10']][[i]]])
}
comuni_sp$M5S_2013_neigh_avg[
  is.na(comuni_sp$M5S_2013_neigh_avg)] <- 0 # Islands

# Correction dist_10k_ppl
camera_sp$dist_10k_ppl <- camera_sp$dist_10k_ppl + 1
comuni_sp$dist_10k_ppl <- comuni_sp$dist_10k_ppl + 1

# Correction income2016_procapite
comuni_sp$income2016_procapite <- 
  comuni_sp$income2016_procapite / 10000

# Make participation binary
camera_sp$meetup2018_5km_buff_freq_7days_bin <- 
  camera_sp$meetup2018_5km_buff_freq_7days>1
comuni_sp$meetup2018_5km_buff_freq_7days_bin <- 
  comuni_sp$meetup2018_5km_buff_freq_7days>1

camera_sp$meetup_treated_2018 <- 
  camera_sp$meetup2018_5km_buff_freq_7days_bin
comuni_sp$meetup_treated_2018 <- 
  comuni_sp$meetup2018_5km_buff_freq_7days_bin

camera_sp$meetup2014_5km_buff_freq_7days_bin <- 
  camera_sp$meetup2014_5km_buff_freq_7days>1
comuni_sp$meetup2014_5km_buff_freq_7days_bin <- 
  comuni_sp$meetup2014_5km_buff_freq_7days>1

camera_sp$meetup_treated_2013 <- 
  camera_sp$meetup2014_5km_buff_freq_7days_bin
comuni_sp$meetup_treated_2013 <- 
  comuni_sp$meetup2014_5km_buff_freq_7days_bin

comuni_sp$meetup2013_5km_buff_freq_7days_bin <- 
  comuni_sp$meetup2014_5km_buff_freq_7days_bin
comuni_sp$meetup2013_5km_buff_freq_7days <- 
  comuni_sp$meetup2014_5km_buff_freq_7days

comuni_sp$meetup2018_treated_diff <- 
  comuni_sp$meetup_treated_2018 + comuni_sp$meetup_treated_2013

comuni_sp$meetup2018_treated_diff <- 
  comuni_sp$meetup_treated_2018 + comuni_sp$meetup_treated_2013
comuni_sp$meetup2018_treated_diff <- 
  comuni_sp$meetup_treated_2018 + comuni_sp$meetup2018_treated_diff

comuni_sp$meetup2018_treated_diff <- 
  factor(comuni_sp$meetup2018_treated_diff, levels = c(0, 1, 3, 2))
levels(comuni_sp$meetup2018_treated_diff) <- 
  c("Never", "2013 only", "2013 & 2018", "2018 only")

sum(comuni_sp$meetup_treated_2018) / nrow(comuni_sp)
sum(comuni_sp$meetup_treated_2013) / nrow(comuni_sp)
sum(comuni_sp$meetup2018_treated_diff == 'Never') / nrow(comuni_sp)
sum(comuni_sp$meetup2018_treated_diff == '2013 & 2018') / nrow(comuni_sp)

comuni_sp$income <- comuni_sp$income2016_procapite

## Add associations
load("secondary_data/istat_censimento_noprofit/01_comuni_harmonized.RData")
comuni_sp <- merge(comuni_sp, comuni_nvol, by = 'PRO_COM_T', all.x = T)
comuni_sp$ANPIV[is.na(comuni_sp$ANPIV)] <- 0 
comuni_sp$ANPINV[is.na(comuni_sp$ANPINV)] <- 0 

comuni_sp$ANPIV_perc <- 
  comuni_sp$ANPIV / comuni_sp$pop2011
```

Note: Turnout for the comune affected by the 2016 earthquake is averaged due to displacements. 

```{r}
## Add referendum turnout
load("secondary_data/minint_referendum_2016/01_comuni2018_turnout2016.RData")
comuni_sp <- 
  merge(comuni_sp, 
        comuni2018_turnout2016[,c('comuni2018',
                                  "registered_voters",
                                  "effective_voters" )], 
        by.x = 'PRO_COM_T', by.y = 'comuni2018', all.x = T)
comuni_sp$ref16_turnout <- 
  comuni_sp$effective_voters / comuni_sp$registered_voters

utm_crs <- 
  '+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs'
lonlat_crs <- 
  '+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0'

earthquake_pnt <- 
  SpatialPoints(cbind(13.22, 42.71),
                proj4string = CRS(lonlat_crs))

buffer_earthquake <- 
  gBuffer(spTransform(earthquake_pnt,
                      utm_crs), width = 50000)
buffer_plot <- 
  gBuffer(spTransform(earthquake_pnt,
                      utm_crs), width = 250000)
res1 <- gWithin(comuni_sp, spTransform(buffer_earthquake, lonlat_crs),
                byid = TRUE)[1,]
res2 <- gWithin(comuni_sp, spTransform(buffer_plot, lonlat_crs),
                byid = TRUE)[1,]

plot(spTransform(buffer_earthquake, lonlat_crs))
plot(comuni_sp[res1,], add=T, col = 'red')
plot(earthquake_pnt, add=T)
plot(comuni_sp[res2,], add=T)

comuni_sp$ref16_turnout[res1] <- 
  sum(comuni_sp$effective_voters[res1]) / sum(comuni_sp$registered_voters[res1])

```

```{r map_comuni_2018_treated_2018_diff, fig.height=12, fig.width=8}
comuni_sf_keep005 <- st_as_sf(comuni_sp_keep005)
comuni_sf_keep005 <- merge(comuni_sf_keep005, comuni_sp, by = 'PRO_COM_T')

ggplot(comuni_sf_keep005) + 
  geom_sf(aes(fill = as.factor(meetup2018_treated_diff)), 
          colour = 'black', size = .02) + 
  geom_sf(data = regione2018_sf, fill = NA, colour = 'black') +
  scale_fill_manual(values = c("white", "#377eb8","#e41a1c","#4daf4a")) + 
  labs(fill='Comune was "treated" by Meetup events' ,
       caption = paste0("number of Comune in 2018: ", nrow(comuni_sp), 
                        "\nnumber of 2018 Comune \"treated\" by events in 2018: ", 
                        sum(comuni_sp$meetup2018_treated), 
                        " (", 
                        round(sum(comuni_sp$meetup2018_treated) / 
                                nrow(comuni_sp) * 100,2), "%)",
                        "\nnumber of 2018 Comune \"treated\" by events in 2013: ", 
                        sum(comuni_sp$meetup2013_treated), 
                        " (", round(sum(comuni_sp$meetup2014_treated) / 
                                      nrow(comuni_sp) * 100,2), "%)",
                        "\nnumber of 2018 Comune \"treated\" by events in both 2013 & 2018: ", 
                        sum(comuni_sp$meetup2018_treated_diff == "2013 & 2018"), 
                        " (", round(sum(comuni_sp$meetup2018_treated_diff == 
                                          "2013 & 2018") / nrow(comuni_sp) * 100,2), 
                        "%)")) +
  theme_bw() + 
  theme(legend.position = c(0.8,0.8))
```

```{r}
comuni_sp_data <- comuni_sp@data
save(comuni_sp_data, file = "replication_file_09_comuni_sp_data.RData")
```

```{r}
load("replication_file_07_province_sp_list.RData")
load("replication_file_09_comuni_sp_data.RData")
```

```{r}
require(ggplot2)
require(dplyr)
require(sf)

prov.sf  <- 
  st_as_sf(province_sp_list[['2012']])

comuni_centroids.sf <-
  comuni_sp %>%
  st_as_sf() %>%
  st_transform(32632) %>%
  st_centroid()

comuni_centroids.sf <- 
  merge(comuni_centroids.sf %>%
          select(PRO_COM_T), 
        comuni_sp_data %>%
          dplyr::select(PRO_COM_T, ANPIV, pop2011, 
                        popdegree2011, pophousewife2011,
                        registered_voters, effective_voters,
                        popworkforce2011, popnotworkforce2011))

comuni_centroids.sf$lat <- 
  st_coordinates(comuni_centroids.sf)[,2]

res <- 
  st_join(comuni_centroids.sf, prov.sf, join = st_within) %>%
  as.data.frame() %>%
  dplyr::filter(!is.na(SIGLA)) %>%
  dplyr::group_by(SIGLA) %>%
  dplyr::summarize(ANPIV_perc = sum(ANPIV) / population[1],
                   degree_perc = sum(popdegree2011) / 
                     sum(pop2011),
                   housewive_perc = sum(pophousewife2011) / 
                     sum(pop2011),
                   ref16_turnout = sum(effective_voters) /
                     sum(registered_voters),
                   partwf_perc = 
                     sum(popworkforce2011) / 
                     sum(popworkforce2011 + popnotworkforce2011),
                   lat = mean(lat))

prov_with_data.sf <- 
  merge(prov.sf, res, by = "SIGLA")

prov_with_data.sf$puntifi1_mean <- 
  region_sp_list[["2013"]]$puntifi1_mean[match(prov_with_data.sf$COD_REG, 
                                               region_sp_list[["2013"]]$COD_REG)]
prov_with_data.sf$puntifi4_mean <- 
  region_sp_list[["2013"]]$puntifi4_mean[match(prov_with_data.sf$COD_REG, 
                                               region_sp_list[["2013"]]$COD_REG)]
prov_with_data.sf$puntifi8_mean <- 
  region_sp_list[["2013"]]$puntifi8_mean[match(prov_with_data.sf$COD_REG, 
                                               region_sp_list[["2013"]]$COD_REG)]
prov_with_data.sf$puntifi9_mean <- 
  region_sp_list[["2013"]]$puntifi9_mean[match(prov_with_data.sf$COD_REG, 
                                               region_sp_list[["2013"]]$COD_REG)]
prov_with_data.sf$puntifi10_mean <- 
  region_sp_list[["2013"]]$puntifi10_mean[match(prov_with_data.sf$COD_REG, 
                                               region_sp_list[["2013"]]$COD_REG)]

load('secondary_data/istisan_reg_2014/01_istisan_reg_2014.RData')
istisan_reg_2014$donor_perc <- 
  istisan_reg_2014$donors / (istisan_reg_2014$pop1865/10000)
istisan_reg_2014$COD_REG <- 
  c(13, 17, 18, 15, 8, 6, 12, 7, 3, 11, 14, 1, 16, 20, 19, 9, 4, 10, 2, 5)

prov_with_data.sf$donor_perc <- 
  istisan_reg_2014$donor_perc[match(prov_with_data.sf$COD_REG, 
                                    istisan_reg_2014$COD_REG)]

prov_with_data.sf <-
  prov_with_data.sf %>%
  dplyr::mutate(pop_density = population / Shape_Area,
                over65_perc = population_over_65 / population,
                foreignpop_perc = foreignpop / population,
                foreignpop_africa_perc = foreignpop_africa / population,
                income = taxable_income / n_taxpayers,
                nat_pol_trust = puntifi4_mean + puntifi1_mean,
                loc_pol_trust = puntifi8_mean + puntifi9_mean + puntifi10_mean)
  
prov_with_data.sf$north <- 
  as.numeric(as.character(prov_with_data.sf$COD_REG)) < 13

prov_with_data.sf$metropolian_city <-
  prov_with_data.sf$SIGLA %in%  
  c("NA", "MI", "RM", "GE", "VE", "CA", "TO", "BA", "CT", "FI", "BO", 
    "PA", "ME", "RC")
```

### Table 1

```{r, results='asis'}
require(lm.beta)
meetup_prov_lm1 <- 
  lm(n_meetup_events_per_10k ~ 
       log(pop_density) + metropolian_city + unemployment + 
       employment + over65_perc + housewive_perc + partwf_perc +
       foreignpop_perc +  foreignpop_africa_perc + income + degree_perc +
       # ref16_turnout + donor_perc + ANPIV_perc + 
       nat_pol_trust + loc_pol_trust,
     data = prov_with_data.sf) 

meetup_prov_lm2 <- 
  lm(n_meetup_events_per_10k ~ 
       log(pop_density) + metropolian_city + unemployment + 
       employment + over65_perc + housewive_perc + partwf_perc +
       foreignpop_perc + foreignpop_africa_perc + 
       income + degree_perc + 
       ref16_turnout + donor_perc + ANPIV_perc +
       nat_pol_trust + loc_pol_trust,
     data = prov_with_data.sf) 

do.call(
  stargazer, 
  c(list(meetup_prov_lm1,
         meetup_prov_lm2,
         coef = list(lm.beta(meetup_prov_lm1)$standardized.coefficients,
                     lm.beta(meetup_prov_lm2)$standardized.coefficients),
         p = list(summary(meetup_prov_lm1)$coefficients[,4],
                  summary(meetup_prov_lm2)$coefficients[,4]),
         type = 'html', out = "table/meetup_prov_lmod.html",
         # type = 'latex', out = "table/meetup_prov_lmod.tex",
         dep.var.labels = "Meetup events per 10k ppl.",
         column.labels = c("Mod. 1 (Beta)", "Mod. 2 (Beta)"),
         omit = c('pop_density', 'metropolian_city',
                  '^employment', 'over65_perc', 'housewive_perc',
                  'foreignpop_perc', 'partwf_perc',
                  'foreignpop_africa_perc', 'loc_pol_trust'),
         order = c('degree_perc', 
                   'unemployment', 'income', 
                   'nat_pol_trust',
                   'ANPIV_perc', 'donor_perc', 'ref16_turnout'),
         covariate.labels = var_dict_labels[
           c('degree_perc', 
             'unemployment', 'income', 
             'nat_pol_trust',
             'ANPIV_perc', 'donor_perc', 'ref16_turnout')]),
    sg_param_intext))

```

### Table 2 (Supplemental online material)

```{r, results='asis'}
do.call(
  stargazer, 
  c(list(meetup_prov_lm1,
         meetup_prov_lm2,
         coef = list(lm.beta(meetup_prov_lm1)$standardized.coefficients,
                     lm.beta(meetup_prov_lm2)$standardized.coefficients),
         p = list(summary(meetup_prov_lm1)$coefficients[,4],
                  summary(meetup_prov_lm2)$coefficients[,4]),
         type = 'latex', out = "table/meetup_prov_lmod_full.tex",
         # type = 'html', out = "table/meetup_prov_lmod_full.html",
         dep.var.labels = "Meetup events per 10k ppl.",
         column.labels = c("Beta", "Beta")
         ,
         covariate.labels = var_dict_labels[
           c("log(pop_density)", "metropolitan_city",
             "unemployment","employment",
             "over65_perc","housewive_perc","partwf_perc",
             "foreignpop_perc",
             "foreignpop_africa_perc", "income", "degree_perc",
             "ref16_turnout", "donor_perc", "ANPIV_perc",
             "nat_pol_trust", "loc_pol_trust")]
         ),
    sg_param_appendix))
```

### Figure 4

```{r meetup_ANPIV_prov_2012, fig.width=9, fig.height=7}
require(scales)
grid.arrange(
  
ggplot(prov_with_data.sf) +
  geom_sf(aes(fill = ANPIV_perc), size = .05) +
  scale_fill_distiller(palette = "Spectral", 
                       labels = function(x){paste0(round(x*100,0),"%")}) +
  theme_bw() +
  theme(axis.title=element_blank(),
        axis.text=element_blank(),
        axis.ticks=element_blank(),
        legend.position = "bottom") +
  labs(fill = "Volunteers"),

ggplot(prov_with_data.sf) +
  geom_sf(aes(fill = n_meetup_events_per_10k), , size = .05) +
  scale_fill_distiller(palette = "Spectral") +
  theme_bw() + 
  theme(axis.title=element_blank(),
        axis.text=element_blank(),
        axis.ticks=element_blank(),
        legend.position = "bottom") +
  labs(fill = "Meetup events per 10k ppl."),
ncol = 2)
```

### Figure 10 (Supplemental online material)

```{r comuni_density_dist_10k, fig.width=9, fig.height=4}
pop_density_df <- 
  data.frame(pop2011 = comuni_sp_data$pop2011,
             pop_density = comuni_sp_data$pop_density,
             dist_10k_ppl = comuni_sp_data$dist_10k_ppl,
             north = comuni_sp_data$north)

pop_density_df$north_south <- 
  ifelse(pop_density_df$north, "North", "South")

grid.arrange(
  ggplot(data=pop_density_df, aes(x=dist_10k_ppl))+
    geom_density(aes(fill=north_south), alpha = .5) + 
    theme(legend.title = element_blank()) +
    labs(x=var_dict_labels['dist_10k_ppl']),
  ggplot(data=pop_density_df, aes(x=log(pop_density)))+
    geom_density(aes(fill=north_south), alpha = .5) + 
    theme(legend.title = element_blank()) +
    labs(x=var_dict_labels['log(pop_density)']))
```


```{r}

## LM 2018
m5s_2018_comuni_lmmod <- 
  lm(data = comuni_sp_data, 
     formula = "M5S_2018~log(pop_density)+log(dist_10k_ppl)+
      unemployment*north+housewive_perc+
      partwf_perc+over65_perc+degree_perc+foreignpop_perc+foreignpop_africa_perc+
      M5S_2018_neigh_avg+meetup_treated_2018+income+ANPIV_perc+ref16_turnout",
     weights = pop2011, na.action=na.exclude)

m5s_2018_comuni_lmmod_nointeraction <- 
  lm(data = comuni_sp_data, 
     formula = "M5S_2018~log(pop_density)+log(dist_10k_ppl)+
      unemployment+north+housewive_perc+
      partwf_perc+over65_perc+degree_perc+foreignpop_perc+foreignpop_africa_perc+
      M5S_2018_neigh_avg+meetup_treated_2018+income+ANPIV_perc+ref16_turnout",
     weights = pop2011, na.action=na.exclude)

m5s_2018_comuni_lmmod_north <- 
  lm(data = comuni_sp_data[comuni_sp_data$north==TRUE,], 
     formula = "M5S_2018~log(pop_density)+log(dist_10k_ppl)+
      unemployment+housewive_perc+
      partwf_perc+over65_perc+degree_perc+foreignpop_perc+foreignpop_africa_perc+
      M5S_2018_neigh_avg+meetup_treated_2018+income+ANPIV_perc+ref16_turnout",
     weights = pop2011, na.action=na.exclude)

m5s_2018_comuni_lmmod_south <- 
  lm(data = comuni_sp_data[!comuni_sp_data$north,], 
     formula = "M5S_2018~log(pop_density)+log(dist_10k_ppl)+
      unemployment+housewive_perc+
      partwf_perc+over65_perc+degree_perc+foreignpop_perc+foreignpop_africa_perc+
      M5S_2018_neigh_avg+meetup_treated_2018+income+ANPIV_perc+ref16_turnout",
     weights = pop2011, na.action=na.exclude)

m5s_diff_comuni_lmmod <- 
  lm(data = comuni_sp_data, 
     formula = "M5S_diff~log(pop_density)+log(dist_10k_ppl)+
      unemployment*north+housewive_perc+
      partwf_perc+over65_perc+degree_perc+foreignpop_perc+foreignpop_africa_perc+
      M5S_2018_neigh_avg+meetup_treated_2018+income+ANPIV_perc+ref16_turnout",
     weights = pop2011,
     na.action = 'na.exclude')

m5s_diff_comuni_lmmod_north <- 
  lm(data = comuni_sp_data[comuni_sp_data$north==TRUE,], 
     formula = "M5S_diff~log(pop_density)+log(dist_10k_ppl)+
      unemployment*north+housewive_perc+
      partwf_perc+over65_perc+degree_perc+foreignpop_perc+foreignpop_africa_perc+
      M5S_2018_neigh_avg+meetup_treated_2018+income+ANPIV_perc+ref16_turnout",
     weights = pop2011,
     na.action = 'na.exclude')

m5s_diff_comuni_lmmod_south <- 
  lm(data = comuni_sp_data[comuni_sp_data$north==FALSE,], 
     formula = "M5S_diff~log(pop_density)+log(dist_10k_ppl)+
      unemployment*north+housewive_perc+
      partwf_perc+over65_perc+degree_perc+foreignpop_perc+foreignpop_africa_perc+
      M5S_2018_neigh_avg+meetup_treated_2018+income+ANPIV_perc+ref16_turnout",
     weights = pop2011,
     na.action = 'na.exclude')

m5s_diff_comuni_lmmod_nointeraction <- 
  lm(data = comuni_sp_data, 
     formula = "M5S_diff~log(pop_density)+log(dist_10k_ppl)+
      unemployment+north+housewive_perc+
      partwf_perc+over65_perc+degree_perc+foreignpop_perc+foreignpop_africa_perc+
      M5S_2018_neigh_avg+meetup_treated_2018+income+ANPIV_perc+ref16_turnout",
     weights = pop2011,
     na.action = 'na.exclude')

comuni_m5s_lmmod_residuals_sf <- 
  comuni_sf_keep005
comuni_m5s_lmmod_residuals_sf$residuals2018 <- 
  resid(m5s_2018_comuni_lmmod)
```

### Table 5 (Supplemental online material)

```{r, results='asis'}
do.call(
  stargazer, 
  c(list(m5s_2018_comuni_lmmod,
         m5s_2018_comuni_lmmod_nointeraction,
         m5s_2018_comuni_lmmod_north,
         m5s_2018_comuni_lmmod_south,
         type = 'latex', out = "table/m5s2018_lmod_full.tex",
         dep.var.labels = "M5S vote 2018, \\%",
         column.labels = c("w/t interaction term (B)", 
                           "w/tout interaction term (B)", 
                           "North only (B)", 
                           "South only (B)"),
         covariate.labels = var_dict_labels[
           c("log(pop_density)","log(dist_10k_ppl)",
             "unemployment","north",
             "housewive_perc","partwf_perc","over65_perc",
             "degree_perc","foreignpop_perc",
             "foreignpop_africa_perc","M5S_2018_neigh_avg",
             "meetup_treated_2018",
             "income","ANPIV_perc",
             "ref16_turnout", "unemployment:north")]),
    sg_param_appendix))
```

```{r, results='asis'}
stargazer(m5s_diff_comuni_lmmod,
          m5s_diff_comuni_lmmod_nointeraction,
          m5s_diff_comuni_lmmod_north,
          m5s_diff_comuni_lmmod_south,
          type = 'latex', out = "table/m5s_diff_lmod_full.tex",
          dep.var.labels = "M5S vote 2018, \\%",
          column.labels = c("w/t interaction term (B)", 
                            "w/tout interaction term (B)", 
                            "North only (B)", 
                            "South only (B)"),
          covariate.labels = var_dict_labels[
            c("log(pop_density)","log(dist_10k_ppl)",
              "unemployment","north",
              "housewive_perc","partwf_perc","over65_perc",
              "degree_perc","foreignpop_perc",
              "foreignpop_africa_perc","M5S_2018_neigh_avg",
              "meetup_treated_2018",
              "income","ANPIV_perc",
              "ref16_turnout", "unemployment:north")],
          column.sep.width = "4pt",
          omit.stat = c('f','ser'),
          report = "vc*",
          single.row = TRUE,
          model.numbers = FALSE,
          model.names = FALSE,
          float = FALSE,
          style = 'ajps')
```

```{r}
## LM 2013

m5s_2013_comuni_lmmod <- 
  lm(data = comuni_sp_data, 
     formula = "M5S_2013~log(pop_density)+log(dist_10k_ppl)+
     unemployment*north+housewive_perc+partwf_perc+over65_perc+
     degree_perc+foreignpop_perc+foreignpop_africa_perc+M5S_2013_neigh_avg+
     meetup_treated_2013+income+ANPIV_perc+ref16_turnout",
     weights = pop2011, na.action = na.exclude)

m5s_2013_comuni_lmmod_nointeraction <- 
  lm(data = comuni_sp_data, 
     formula = "M5S_2013~log(pop_density)+log(dist_10k_ppl)+
     unemployment+north+housewive_perc+partwf_perc+over65_perc+
     degree_perc+foreignpop_perc+foreignpop_africa_perc+M5S_2013_neigh_avg+
     meetup_treated_2013+income+ANPIV_perc+ref16_turnout",
     weights = pop2011, na.action = na.exclude)

m5s_2013_comuni_lmmod_north <- 
  lm(data = comuni_sp_data[comuni_sp_data$north == TRUE,], 
     formula = "M5S_2013~log(pop_density)+log(dist_10k_ppl)+
     unemployment+housewive_perc+partwf_perc+over65_perc+
     degree_perc+foreignpop_perc+foreignpop_africa_perc+M5S_2013_neigh_avg+
     meetup_treated_2013+income+ANPIV_perc+ref16_turnout",
     weights = pop2011, na.action = na.exclude)

m5s_2013_comuni_lmmod_south <- 
  lm(data = comuni_sp_data[comuni_sp_data$north == FALSE,], 
     formula = "M5S_2013~log(pop_density)+log(dist_10k_ppl)+
     unemployment+housewive_perc+partwf_perc+over65_perc+
     degree_perc+foreignpop_perc+foreignpop_africa_perc+M5S_2013_neigh_avg+
     meetup_treated_2013+income+ANPIV_perc+ref16_turnout",
     weights = pop2011, na.action = na.exclude)

comuni_m5s_lmmod_residuals_sf$residuals2013 <- 
  resid(m5s_2013_comuni_lmmod)
```

### Table 3

```{r, results = 'asis'}
do.call(
  stargazer, c(list(m5s_2013_comuni_lmmod,
                    m5s_2018_comuni_lmmod,
                    # type = 'latex', out = "table/m5s2018-13_lmod.tex",
                    type = 'html', out = "table/m5s2018-13_lmod.html",
                    dep.var.labels = "M5S vote, \\%", 
                    column.labels = c("2013 (B)", "2018 (B)"),
                    omit = c('^northTRUE$', 'partwf_perc', 'pop\\_density', 'over65', 
                             'foreignpop', 'degree', 'year',
                             'Constant', "dist\\_10k_ppl","housewive"),
                    order = c("^unemployment$",
                              "meetup_treated_2013", "meetup_treated_2018",
                              "income", "ANPIV_perc", "ref16_turnout",
                              "M5S\\_2013\\_neigh", "M5S\\_2018\\_neigh", "unemployment:north"),
                    covariate.labels = c("Unemployment (pct) | North = 0", var_dict_labels[
                      c("meetup_treated_2013",
                        "meetup_treated_2018",
                        "income","ANPIV_perc","ref16_turnout", 
                        "M5S_2013_neigh_avg", "M5S_2018_neigh_avg", "unemployment:north")
                      ])),
               sg_param_intext))
```

### Table 4 (Supplemental online material)

```{r, results = 'asis'}
do.call(
  stargazer, 
  c(list(m5s_2013_comuni_lmmod,
         m5s_2013_comuni_lmmod_nointeraction,
         m5s_2013_comuni_lmmod_north,
         m5s_2013_comuni_lmmod_south,
         type = 'latex', out = "table/m5s2013_lmod_full.tex",
         dep.var.labels = "M5S vote 2013, \\%",
         column.labels = c("w/t interaction term (B)", 
                           "w/tout interaction term (B)", 
                           "North only (B)", 
                           "South only (B)"),
         covariate.labels = var_dict_labels[
           c("log(pop_density)","log(dist_10k_ppl)",
             "unemployment","north",
             "housewive_perc","partwf_perc","over65_perc",
             "degree_perc","foreignpop_perc",
             "foreignpop_africa_perc","M5S_2013_neigh_avg",
             "meetup_treated_2013",
             "income","ANPIV_perc",
             "ref16_turnout", "unemployment:north")
           ]),
    sg_param_appendix))
```

### Figure 9 (Supplemental online material)

```{r m5s2018-13_lmod_residuals, fig.width=9, fig.height=6}
grid.arrange(
    ggplot(comuni_m5s_lmmod_residuals_sf, aes(fill = residuals2013)) + 
    geom_sf(colour = NA) + scale_fill_viridis() +
    theme_bw() + theme(legend.position = 'bottom') + labs(fill = '2013'),
  ggplot(comuni_m5s_lmmod_residuals_sf, aes(fill = residuals2018)) + 
    geom_sf(colour = NA) + scale_fill_viridis() +
    theme_bw() + theme(legend.position = 'bottom') + labs(fill = '2018'),
  ncol = 2)
```


### Figure 3

```{r meetup_events_freq_2018-13, fig.height=6, fig.width=6}
grid.arrange(
  ggplot(data=comuni_sp_data, aes(x=as.factor(meetup2018_5km_buff_freq_7days))) + 
    geom_bar(aes(y = (..count..)/sum(..count..))) + 
    theme_bw() + scale_y_continuous(labels = percent) + 
    labs(x = "6 Jan - 3 Mar 2018", y = NULL) +
    geom_segment(aes(x = 3, xend = 10, y = .2, yend = .2), linetype = 2) +
    geom_text(aes(y = .25, x = 6, label = 'Treatement')),
  ggplot(data=comuni_sp_data, aes(x=as.factor(meetup2014_5km_buff_freq_7days))) + 
    geom_bar(aes(y = (..count..)/sum(..count..))) + 
    theme_bw() + scale_y_continuous(labels = percent) + 
    labs(x = "6 Jan - 24 Feb 2013", y = NULL) +
    geom_segment(aes(x = 3, xend = 9, y = .2, yend = .2), linetype = 2) +
    geom_text(aes(y = .25, x = 6, label = 'Treatement')),
  ncol = 1)
```

```{r, results = 'asis'}
load('replication_file_10_camera_2013_by_region.RData')

camera_2013_by_region_df$DEN_REG <- 
  region_sp_list$`2013`$DEN_REG[
    match(camera_2013_by_region_df$COD_REG, region_sp_list$`2013`$COD_REG)]

reg_long_df <- data.frame()

for (y in 2013:2016) {
  reg_long_df <- 
    rbind(reg_long_df,
          data.frame(name = 
                       region_sp_list[[as.character(y)]]$DEN_REG,
                     COD_REG =
                       region_sp_list[[as.character(y)]]$COD_REG,
                     n_meetup_events = 
                       region_sp_list[[as.character(y)]]$n_meetup_events,
                     population = 
                       region_sp_list[[as.character(y)]]$population,
                     population_over_65 = 
                       region_sp_list[[as.character(y)]]$population_over_65,
                     unemployment = 
                       region_sp_list[[as.character(y)]]$unemployment,
                     employment = 
                       region_sp_list[[as.character(y)]]$employment,
                     foreignpop = 
                       region_sp_list[[as.character(y)]]$foreignpop,
                     foreignpop_africa = 
                       region_sp_list[[as.character(y)]]$foreignpop_africa,
                     income = 
                       region_sp_list[[as.character(y)]]$income/10000,
                     
                     ppapo_mean = 
                       region_sp_list[[as.character(y)]]$ppapo_mean,
                     psind_mean = 
                       region_sp_list[[as.character(y)]]$psind_mean,
                     pgrvo_mean = 
                       region_sp_list[[as.character(y)]]$pgrvo_mean,
                     paeco_mean = 
                       region_sp_list[[as.character(y)]]$paeco_mean,
                     
                     puntifi1_mean = 
                       region_sp_list[[as.character(y)]]$puntifi1_mean,
                     puntifi2_mean = 
                       region_sp_list[[as.character(y)]]$puntifi2_mean,
                     puntifi3_mean = 
                       region_sp_list[[as.character(y)]]$puntifi3_mean,
                     puntifi4_mean = 
                       region_sp_list[[as.character(y)]]$puntifi4_mean,
                     puntifi5_mean = 
                       region_sp_list[[as.character(y)]]$puntifi5_mean,
                     puntifi8_mean = 
                       region_sp_list[[as.character(y)]]$puntifi8_mean,
                     puntifi9_mean = 
                       region_sp_list[[as.character(y)]]$puntifi9_mean,
                     puntifi10_mean = 
                       region_sp_list[[as.character(y)]]$puntifi10_mean,
                     
                     pol_trust_diff_mean =
                       region_sp_list[[
                         as.character(y)]]$pol_trust_diff_mean,
                     
                     comiz_mean = 
                       region_sp_list[[as.character(y)]]$comiz_mean,
                     cortei_mean = 
                       region_sp_list[[as.character(y)]]$cortei_mean,
                     dibpo_mean = 
                       region_sp_list[[as.character(y)]]$dibpo_mean,
                     
                     volon_mean = 
                       region_sp_list[[as.character(y)]]$volon_mean,
                     degree_perc = 
                       region_sp_list[[as.character(y)]]$degree_perc,
                     
                     area = 
                       region_sp_list[[as.character(y)]]$Shape_Area,
                     year = y))
}

## Add blood donation
load('secondary_data/istisan_reg_2014/01_istisan_reg_2014.RData')
istisan_reg_2014$donor_perc <- 
  istisan_reg_2014$donors / (istisan_reg_2014$pop1865/10000)

reg_long_df <- 
  merge(reg_long_df, 
        istisan_reg_2014[,c('Regione','donor_perc')], by.x='name', by.y="Regione")

## Add associations
load('secondary_data/istat_censimento_noprofit/01_comuni_harmonized.RData')

regioni_nvol <-
  comuni_nvol %>%
  dplyr::group_by(DEN_REG) %>%
  dplyr::summarize(ANPIV = sum(ANPIV, na.rm = T))

reg_long_df <- 
  merge(reg_long_df, regioni_nvol[,c('DEN_REG', "ANPIV")],
        by.x='name', by.y="DEN_REG")

reg_long_df$ANPIV_perc <- 
  reg_long_df$ANPIV / reg_long_df$population

## Add referendum turnout
load("secondary_data/minint_referendum_2016/01_comuni2018_turnout2016.RData")
regioni_turnout <- 
  comuni2018_turnout2016 %>%
  dplyr::group_by(region_nam) %>%
  dplyr::summarize(ref16_turnout = 
                     sum(effective_voters) / sum(registered_voters))

reg_long_df <- 
  merge(reg_long_df, 
        regioni_turnout[,c('region_nam', "ref16_turnout")],
        by.x='name', by.y="region_nam")

assignMacroRegion <- function(x) {
  if (x %in% c('Piemonte','Lombardia',
               'Liguria',"Valle d'Aosta/Vallée d'Aoste")) {
    return('North-west')
  }
  if(x %in% c('Veneto','Friuli-Venezia Giulia',
              "Trentino-Alto Adige/Südtirol",
              "Emilia-Romagna")) {
    return('North-east')
  }
  if (x %in% c("Toscana", "Umbria", "Marche", "Lazio")) {
    return('Centre')
  }
  if (x %in% c("Abruzzo", "Puglia", "Basilicata", 
               "Campania", "Calabria", "Molise")) {
    return('South')
  }
  if (x %in% c("Sicilia", "Sardegna")) {
    return("Islands")
  }
}

reg_long_df$macro_region <- 
  sapply(reg_long_df$name, assignMacroRegion)
reg_long_df$north_south <- 
  ifelse(reg_long_df$macro_region %in% 
           c('South','Islands'), "South", "North")
reg_long_df$north <- 
  reg_long_df$north_south == "North"
reg_long_df$n_meetup_events_per_10k <- 
  reg_long_df$n_meetup_events / (reg_long_df$population/10000)

reg_long_df <- 
  reg_long_df %>%
  dplyr::group_by(name) %>%
  dplyr::arrange(year, .by_group = TRUE) %>%
  dplyr::mutate(pop_density = population / area,
                over65_perc = population_over_65 / population,
                foreignpop_perc = foreignpop / population,
                foreignpop_africa_perc = foreignpop_africa / population)

reg_long_df$nat_pol_trust <- 
  with(reg_long_df, puntifi4_mean + puntifi1_mean)
reg_long_df$loc_pol_trust <- 
  with(reg_long_df, puntifi8_mean + puntifi9_mean + puntifi10_mean)
reg_long_df$pol_trust_diff <- 
  reg_long_df$pol_trust_diff_mean

reg_long_df$veneto <- 
  as.numeric(reg_long_df$name == 'Veneto')

reg_2013_long_df <- 
  reg_long_df[reg_long_df$year == 2013,]
reg_2013_long_df$M5S_2013 <-
  camera_2013_by_region_df$M5S_2013[
    match(reg_2013_long_df$name, 
          camera_2013_by_region_df$DEN_REG)]

reg_long_df$unemployment <- reg_long_df$unemployment / 100
reg_long_df$employment <- reg_long_df$employment / 100

meetup_reg_lmmod1 <-  
  lm(n_meetup_events_per_10k ~ 
       log(pop_density) + unemployment * north + 
       employment + over65_perc + foreignpop_perc + 
       foreignpop_africa_perc + income + degree_perc + 
       nat_pol_trust + loc_pol_trust + 
       year, 
     data = reg_long_df)

meetup_reg_lmmod2 <-  
  lm(n_meetup_events_per_10k ~ log(pop_density) + unemployment * north + 
       employment + over65_perc + foreignpop_perc + 
       foreignpop_africa_perc + income + degree_perc + 
       nat_pol_trust + loc_pol_trust + 
       donor_perc + ANPIV_perc + ref16_turnout + 
       year, data = reg_long_df)

# Predict with interaction effect
newdata <- 
  reg_long_df %>%
          dplyr::group_by(north) %>%
          dplyr::select(north, pop_density, 
                        employment, over65_perc, foreignpop_perc, 
                        foreignpop_africa_perc, income, degree_perc, 
                        nat_pol_trust, loc_pol_trust, 
                        donor_perc, ANPIV_perc, ref16_turnout, year) %>%
          dplyr::summarize_all(median)

newdata <- 
  newdata[c(rep(1, 100),
            rep(2, 100)), ]
newdata$unemployment <- 
  c(seq(from = min(reg_long_df$unemployment[reg_long_df$north == FALSE]),
      to  = max(reg_long_df$unemployment[reg_long_df$north == FALSE]),
      length.out = 100),
    seq(from = min(reg_long_df$unemployment[reg_long_df$north == TRUE]),
      to  = max(reg_long_df$unemployment[reg_long_df$north == TRUE]),
      length.out = 100))

predict.df <- 
  predict(meetup_reg_lmmod2, 
        newdata, 
        type="response", interval="confidence") %>%
  as.data.frame()
predict.df$unemployment <-
  newdata$unemployment
predict.df$north <-
  newdata$north

# Max difference
max_unemp_north <-  
  with(predict.df, max(unemployment[north == TRUE]))

predict_north <-
  predict.df$fit[predict.df$unemployment == max_unemp_north]

predict_south <-
  predict.df$fit[with(predict.df,
                    which.min(abs(unemployment[north == FALSE] - 
                                    max_unemp_north)))]
# Difference north/south when unemployment is about:
max_unemp_north

predict_north - predict_south

meetup_reg_femod1 <- 
  plm::plm(n_meetup_events_per_10k ~ log(pop_density) + 
        unemployment * north + employment + over65_perc + foreignpop_perc + 
        foreignpop_africa_perc + income + degree_perc + 
        nat_pol_trust + loc_pol_trust,
      index = c("year","name"), model = "within", 
      data = as.data.frame(reg_long_df))

meetup_reg_femod2 <- 
  plm(n_meetup_events_per_10k ~ log(pop_density) + unemployment * north + 
        employment + over65_perc + foreignpop_perc + 
        foreignpop_africa_perc + income + degree_perc + 
        nat_pol_trust + loc_pol_trust + 
        donor_perc + ANPIV_perc + ref16_turnout,
      index = c("year","name"), model = "within", 
      data = as.data.frame(reg_long_df))

meetup_reg_lmmod1_stdcoef <- lm.beta(meetup_reg_lmmod1)
meetup_reg_lmmod2_stdcoef <- lm.beta(meetup_reg_lmmod2)

do.call(stargazer,
        c(list(meetup_reg_lmmod1,
               meetup_reg_lmmod2,
               type = 'latex', out = "table/meetup_reg_lmod.tex",
               # type = 'html', out = "table/meetup_reg_lmod.html",
               coef = list(meetup_reg_lmmod1_stdcoef$standardized.coefficients,
                           meetup_reg_lmmod2_stdcoef$standardized.coefficients),
               p = list(summary(meetup_reg_lmmod1)$coefficients[,4],
                        summary(meetup_reg_lmmod2)$coefficients[,4]),
               dep.var.labels = "Meetup events per 10k ppl.",
               column.labels = c("Mod. 1 (Beta)", "Mod. 2 (Beta)"),
               omit = c('density', 'over65', 'foreignpop', 'year', 'Constant', 'degree', '^north'),
               order=c("^unemployment$",
                       "^employment$", "income",
                       "^nat_pol","^loc_pol",
                       "donor_perc", "ANPIV_perc", "ref16_turnout",
                       "unemployment:north"),
               covariate.labels = c("Unemployment (pct) | North = 0", var_dict_labels[
                 c('employment','income',
                   "nat_pol_trust","loc_pol_trust",
                   "donor_perc","ANPIV_perc", "ref16_turnout",
                   "unemployment:north")])),
          sg_param_intext))
  
```

```{r predict_reg_interaction, fig.width=6, fig.height=6}
ggplot() +
  geom_line(data = predict.df, 
            aes(y=fit, 
                x = unemployment,
                colour = north)) +
  geom_ribbon(data = predict.df, 
              aes(x=unemployment, ymax=upr, ymin=lwr,fill=north),
              alpha = 0.2) + 
  geom_point(data =
               reg_long_df, aes(colour=north, x=unemployment,
                                y=n_meetup_events_per_10k,shape=north)) +
  theme_bw() + scale_x_continuous(label = percent) +
  scale_color_manual(values = two_colours) +
  scale_fill_manual(values = two_colours) +
  scale_shape_discrete(labels = c("Southern region", "Northern region")) +
  theme(legend.position = c(0.7,0.7)) +
  labs(x = var_dict_labels['unemployment'], y = "Meetup events per 10k ppl.", shape = NULL) + 
  guides(fill = FALSE, colour = FALSE) +
  geom_text(data = data.frame(x=c(max(predict.df$unemployment[predict.df$north==TRUE])+.01,
                                  max(predict.df$unemployment[predict.df$north==FALSE])+.01),
                              y=c(max(predict.df$fit[predict.df$north==TRUE]),
                                  min(predict.df$fit[predict.df$north==FALSE])),
                              text=c("North","South")),
            aes(x=x,y=y,label=text))
```


```{r meetup_reg_lmmod_residuals, fig.width=10, fig.height=3.75}

# residuals plot
meetup_reg_lmmod_residuals <- 
  data.frame(residuals = summary(meetup_reg_lmmod2)$residuals,
             year = reg_long_df$year,
             region = reg_long_df$name)

reg_residuals_2013_sf <- 
  merge(regione2018_sf,
        meetup_reg_lmmod_residuals[meetup_reg_lmmod_residuals$year == 2013,],
        by.x = 'DEN_REG', by.y = 'region')
reg_residuals_2014_sf <- 
  merge(regione2018_sf, 
        meetup_reg_lmmod_residuals[meetup_reg_lmmod_residuals$year == 2014,],
        by.x = 'DEN_REG', by.y = 'region')
reg_residuals_2015_sf <- 
  merge(regione2018_sf, 
        meetup_reg_lmmod_residuals[meetup_reg_lmmod_residuals$year == 2015,],
        by.x = 'DEN_REG', by.y = 'region')
reg_residuals_2016_sf <- 
  merge(regione2018_sf, 
        meetup_reg_lmmod_residuals[meetup_reg_lmmod_residuals$year == 2016,],
        by.x = 'DEN_REG', by.y = 'region')

grid.arrange(
  ggplot(reg_residuals_2013_sf, aes(fill = residuals)) + 
    geom_sf() + scale_fill_viridis() + theme_bw() + theme(legend.position = 'bottom') +
    labs(fill = "2013") +
    theme(axis.title=element_blank(),
          axis.text=element_blank(),
          axis.ticks=element_blank()),
  ggplot(reg_residuals_2014_sf, aes(fill = residuals)) + 
    geom_sf() + scale_fill_viridis() + theme_bw() + theme(legend.position = 'bottom') +
    labs(fill = "2014") +
    theme(axis.title=element_blank(),
          axis.text=element_blank(),
          axis.ticks=element_blank()),
  ggplot(reg_residuals_2015_sf, aes(fill = residuals)) + 
    geom_sf() + scale_fill_viridis() + theme_bw() + theme(legend.position = 'bottom') +
    labs(fill = "2015") +
    theme(axis.title=element_blank(),
          axis.text=element_blank(),
          axis.ticks=element_blank()),
  ggplot(reg_residuals_2016_sf, aes(fill = residuals)) + 
    geom_sf() + scale_fill_viridis() + theme_bw() + theme(legend.position = 'bottom') +
    labs(fill = "2016") +
    theme(axis.title=element_blank(),
          axis.text=element_blank(),
          axis.ticks=element_blank()),
  ncol = 4)
```

```{r, results = 'asis'}
do.call(
  stargazer, c(
    list(meetup_reg_lmmod1,
          meetup_reg_lmmod1,
          meetup_reg_lmmod2,
          meetup_reg_lmmod2,
          meetup_reg_femod1,
          meetup_reg_femod2,
          type = 'latex', out = "table/meetup_reg_lmod_femod_full.tex",
          coef = list(meetup_reg_lmmod1$coefficients,
                      meetup_reg_lmmod1_stdcoef$standardized.coefficients,
                      meetup_reg_lmmod2$coefficients,
                      meetup_reg_lmmod2_stdcoef$standardized.coefficients,
                      meetup_reg_femod1$coefficients,
                      meetup_reg_femod2$coefficients),
          p = list(summary(meetup_reg_lmmod1)$coefficients[,4],
                   summary(meetup_reg_lmmod1)$coefficients[,4],
                   summary(meetup_reg_lmmod2)$coefficients[,4],
                   summary(meetup_reg_lmmod2)$coefficients[,4],
                   summary(meetup_reg_femod1)$coefficients[,4],
                   summary(meetup_reg_femod2)$coefficients[,4]),
          dep.var.labels = "Meetup events (per 10k people)",
          column.labels = c("(B)", "(beta)", 
                            "(B)", "(beta)",
                            "(B)", "(B)"),
          order = c('north', 'log(pop_density)',
                    'unemployment', '^employment', 
                    'over65_perc', 
                    'foreignpop_perc',
                'foreignpop_africa_perc',
                'income', 'degree_perc', 'nat_pol_trust',
                'loc_pol_trust', 'donor_perc',
                'ANPIV_perc', 'ref16_turnout',
                'year', 'unemployment:north'),
                    covariate.labels =
            c("North (bin) | Unemployment = 0",
              var_dict_labels['log(pop_density)'],
              "Unemployment (pct) | North = 0",
              var_dict_labels[
              c('^employment', 
                    'over65_perc', 
                    'foreignpop_perc',
                'foreignpop_africa_perc',
                'income', 'degree_perc', 'nat_pol_trust',
                'loc_pol_trust', 'donor_perc',
                'ANPIV_perc', 'ref16_turnout',
                'year', 'unemployment:north')])),
         sg_param_appendix))
```

## Matching

```{r}
var <- 
  c('PRO_COM_T', 'M5S_2018','M5S_2013','pop_density', 
    'dist_10k_ppl', 'degree_perc', 'meetup_treated_2018', 
    'meetup_treated_2013', 'unemployment', 'housewive_perc', 
    'partwf_perc', 'over65_perc', "foreignpop_perc", 
    "foreignpop_africa_perc", "income2016_procapite", "M5S_2018_neigh_avg",
    "M5S_2013_neigh_avg", "pop2011",
    "macro_region", "lon", "lat", "M5S_diff", "north", 
    "ANPIV_perc", "ref16_turnout")

comuni_df <- 
  comuni_sp_data %>%  
  dplyr::select(one_of(var)) %>%
  na.omit()

comuni_df$income <- 
  comuni_df$income2016_procapite

# Recode
comuni_df$income2016_procapite_quintile <- 
  cut(comuni_df$income2016_procapite, 
      breaks = quantile(comuni_df$income2016_procapite, 
                        probs = seq(0,1, by=.2)),
      include.lowest = T)


## 2013
m.out <- 
  matchit(meetup_treated_2013 ~ 
            income2016_procapite_quintile + dist_10k_ppl, 
          data = comuni_df, 
          method = "exact")
m.data <- match.data(m.out)

match_2013_residuals_sf <- 
  comuni_sf_keep005
match_2013_residuals_sf <- 
  match_2013_residuals_sf[match_2013_residuals_sf$PRO_COM_T %in% 
                            m.data$PRO_COM_T,]

m5s_2013_comuni_match_lmmod <-
  lm(M5S_2013 ~ meetup_treated_2013 + log(pop_density) + 
       unemployment * north + housewive_perc + 
       partwf_perc + over65_perc + degree_perc + 
       foreignpop_perc + foreignpop_africa_perc + income + 
       M5S_2013_neigh_avg + ANPIV_perc + ref16_turnout,
     weights = pop2011,
     data = m.data, na.action = na.exclude)

match_2013_residuals_sf$residuals <- 
  resid(m5s_2013_comuni_match_lmmod)
match_2013_residuals_sf$meetup_treated_2013 <-
  m.data$meetup_treated_2013

## 2018
m.out <- 
  matchit(meetup_treated_2018 ~ 
            income2016_procapite_quintile + dist_10k_ppl, 
          data = comuni_df, 
          method = "exact")
m.data <- match.data(m.out)

m5s_2018_comuni_match_lmmod <-
  lm(M5S_2018 ~ meetup_treated_2018 + log(pop_density) + unemployment * north  + 
       housewive_perc + partwf_perc + over65_perc + degree_perc + foreignpop_perc +
       foreignpop_africa_perc + income + M5S_2018_neigh_avg + ANPIV_perc + 
       ref16_turnout,
     weights = pop2011,
     data = m.data, na.action = na.exclude)

match_2018_residuals_sf <- comuni_sf_keep005
match_2018_residuals_sf <- 
  match_2018_residuals_sf[match_2018_residuals_sf$PRO_COM_T %in% 
                            m.data$PRO_COM_T,]

match_2018_residuals_sf$residuals <- 
  resid(m5s_2018_comuni_match_lmmod)

match_2018_residuals_sf$meetup_treated_2018 <-
  m.data$meetup_treated_2018

m5s_diff_comuni_match_lmmod <-
  lm(M5S_diff ~ meetup_treated_2018 + log(pop_density) + unemployment * north  + 
       housewive_perc + partwf_perc + over65_perc + degree_perc + 
       foreignpop_perc + foreignpop_africa_perc + income + M5S_2018_neigh_avg  + 
       ANPIV_perc + ref16_turnout,
     weights = pop2011,
     data = m.data)
```

```{r, results = 'asis'}
stargazer(m5s_2013_comuni_match_lmmod,
          m5s_2018_comuni_match_lmmod,
          type = 'latex', out = "table/m5s_matching_lmod.tex",
          dep.var.labels = "M5S vote",
          column.labels = c("\\%, 2013", "\\%, 2018"),
          omit = c('pop\\_density', 'over65', 'foreignpop', 'degree', 'year', 
                   'Constant', "dist\\_10k_ppl","housewive","partwf",
                   "M5S\\_2018\\_neigh", "M5S\\_2013\\_neigh"),
          omit.stat = c('adj.rsq', 'f', 'ser'),
          covariate.labels = 
            var_dict_labels[
              c('meetup_treated_2013', 'meetup_treated_2018',
                'unemployment', 'north', 'income',
                'ANPIV_perc', 'ref16_turnout', 
                'unemployment:north')],
          report = "vc*",
          model.numbers = FALSE,
          model.names = FALSE,
          float = FALSE,
          style = 'ajps')
```

```{r, results = 'asis'}
stargazer(m5s_2013_comuni_match_lmmod,
          m5s_2018_comuni_match_lmmod,
          m5s_diff_comuni_match_lmmod,
          type = 'latex', out = "table/m5s_matching_lmod_full.tex",
          dep.var.labels = "M5S vote",
          column.labels = c("\\%, 2013", "\\%, 2018", "diff. 2018/13"),
          covariate.labels = var_dict_labels[
            c('meetup_treated_2013', 'meetup_treated_2018',
              'log(pop_density)', 'unemployment', 'north',
              'housewive_perc', 'partwf_perc', 'over65_perc',
              'degree_perc', 'foreignpop_perc', 'foreignpop_africa_perc',
              'income', 'M5S_2013_neigh_avg', 'M5S_2018_neigh_avg',
              'ANPIV_perc', 'ref16_turnout', 'unemployment:north')
            ],
          report = "vc*",
          model.numbers = FALSE,
          model.names = FALSE,
          float = FALSE,
          style = 'ajps')
```

```{r comuni_matched_treated_residuals_2013, fig.width=9, height=4.8}
grid.arrange(
  ggplot(match_2013_residuals_sf, aes(fill = as.factor(meetup_treated_2013))) +
    scale_fill_manual(values = c('#66c2a5','#fc8d62'), name = 'Treated in 2013', 
                      labels = c("No", "Yes")) +
    geom_sf(size = .05, colour = 'black') + theme_bw() + 
    geom_sf(data = regione2018_sf, fill = NA) +
    theme(legend.position = 'bottom') +
    labs(caption = paste0('Matched comune (n=', nrow(match_2013_residuals_sf),', 2013)')),
  ggplot(match_2013_residuals_sf, aes(fill = residuals)) +
    scale_fill_viridis(name = 'Residuals from OLS model') +
    geom_sf(size = .05, colour = 'black') + theme_bw() + 
    geom_sf(data = regione2018_sf, fill = NA) +
    theme(legend.position = 'bottom'),
  ncol = 2)
```

```{r comuni_matched_treated_residuals_2018, width=9, height=4.8}
grid.arrange(
  ggplot(match_2018_residuals_sf, aes(fill = as.factor(meetup_treated_2018))) +
    scale_fill_manual(values = c('#66c2a5','#fc8d62'),
                      name = 'Treated in 2018', labels = c("No", "Yes")) +
    geom_sf(size = .05, colour = 'black') + theme_bw() + 
    geom_sf(data = regione2018_sf, fill = NA) +
    theme(legend.position = 'bottom') +
    labs(caption = paste0('Matched comune (n=', 
                          nrow(match_2018_residuals_sf),', 2018)')),
  ggplot(match_2018_residuals_sf, aes(fill = residuals)) +
    scale_fill_viridis(name = 'Residuals from OLS model') +
    geom_sf(size = .05, colour = 'black') + theme_bw() + 
    geom_sf(data = regione2018_sf, fill = NA) +
    theme(legend.position = 'bottom'),
  ncol = 2)
```

## Individual level models 

```{r, results = 'asis'}
itanes_2013 <- 
  read_spss("secondary_data/itanes/ITA2013_(envers2015_01_19).sav")
itanes_2018 <- 
  read_spss("secondary_data/itanes/2019_03_27_Itanes Pastel 2018_release_01_pre-post.sav")

itanes_2013$M5S <-  as.integer(itanes_2013$d90==5)
itanes_2018$M5S <- as.integer(itanes_2018$votoC1==4)

# Closeness to M5S
recodeD43 <- function(x) {
  x <- as.numeric(x)
  x[x %in% c(4,5)] <- NA
  return(x)
}
itanes_2013$clsness_M5S <- 
  ifelse(itanes_2013$d41 == 12, 1, 0) 
itanes_2013$clsness_M5S <-
  itanes_2013$clsness_M5S * recodeD43(itanes_2013$d43)

itanes_2013$unemployed <- itanes_2013$d4 %in% c(5,6,8)
itanes_2018$unemployed <- itanes_2018$d10_7 %in% c(4:6,9)
itanes_2013$degree <- itanes_2013$stu %in% 9:10
itanes_2018$degree <- itanes_2018$d10_11 %in% 8:11
itanes_2013$north <- itanes_2013$regint < 13
itanes_2018$north <- itanes_2018$AreaQ < 4

# Economic stress
itanes_2013$eco_stress <- 
  rescale(
  1*(itanes_2013$lav_1 == 1) + 1*(itanes_2013$lav_2 == 1) +
  0.5*(itanes_2013$lav2 == 1) + 0.25*(itanes_2013$lav2 == 2) +
  0.75*(itanes_2013$lav2b1 == 1) + 0.75*(itanes_2013$lav2b2 == 1) +
  1*(itanes_2013$lav3 == 1) + 0.5*(itanes_2013$lav3 == 2),
  to = c(0,10))

itanes_2013$eco_perception <- 
  rescale(rowMeans(cbind(
    recode(as.numeric(as.character(itanes_2013$d5)), 
           '1' = 1, 
           '2' = 2, 
           '3' = 3,
           '4' = 4,
           '5' = 5,
           '6' = 3)*1,
    recode(as.numeric(as.character(itanes_2013$d5)), 
           '1' = 1, 
           '2' = 2, 
           '3' = 3,
           '4' = 4,
           '5' = 5,
           '6' = 3)*.75
  )), to = c(10,0))

itanes_2018$eco_perception <- 
  rescale(rowMeans(cbind(
    recode(as.numeric(as.character(itanes_2018$d1_1)), 
           '1' = 1, 
           '2' = 2, 
           '3' = 3,
           '4' = 4,
           '5' = 5,
           '6' = 3)*1,
    recode(as.numeric(as.character(itanes_2018$d1_2)), 
           '1' = 1, 
           '2' = 2, 
           '3' = 3,
           '4' = 4,
           '5' = 5,
           '6' = 3)*.75
  )), to = c(10,0))

# Political partecipation
itanes_2013$pol_part <- 
  rescale(
    with(itanes_2013,
         (d78_8 == 1) +
           (d21_1 == 1) + 
           (d21_2 == 1) +
           (d21_3 == 1) +
           (d21_4 == 1) +
           (d21_5 == 1) +
           (d21_6 == 1) +
           (d21_7 == 1) +
           (d21_8 == 1)),
    to = c(0,10))
  

# Online political discussion
itanes_2013$online_pol_discussion <- 0
itanes_2013$online_pol_discussion[itanes_2013$d72_5 == 3] <- 1
itanes_2013$online_pol_discussion[itanes_2013$d72_5 == 2] <- 2
itanes_2013$online_pol_discussion[itanes_2013$d72_5 == 1] <- 3
itanes_2013$online_pol_discussion[itanes_2013$d72_5 == 6] <- NA

# INFORMATION
itanes_2013$info_internet <- 
  itanes_2013$d57 == 7 | itanes_2013$d57b == 7

itanes_2018$info_internet_news <- 
  itanes_2018$d4_1 == 4 | itanes_2018$d4_1_1 == 4
itanes_2018$info_social_media <- 
  itanes_2018$d4_1 == 5 | itanes_2018$d4_1_1 == 5

itanes_2013$info_personal <- 
  itanes_2013$d57 == 2 | itanes_2013$d57b == 2 |
  itanes_2013$d57 == 6 | itanes_2013$d57b == 6

itanes_2018$info_personal <- 
  itanes_2018$d4_1 == 6 | itanes_2018$d4_1_1 == 6

itanes_2013$info_written <- 
  itanes_2013$d57 == 3 | itanes_2013$d57b == 3 |
  itanes_2013$d57 == 4 | itanes_2013$d57b == 4

itanes_2018$info_written <- 
  itanes_2018$d4_1 == 2 | itanes_2018$d4_1_1 == 2

itanes_2013$info_tv <- 
  itanes_2013$d57 == 1 | itanes_2013$d57b == 1 

itanes_2018$info_tv <- 
  itanes_2018$d4_1 == 1 | itanes_2018$d4_1_1 == 1

itanes_2013$internet_mediated_part <- 
  as_factor(itanes_2013$d72_6)
itanes_2013$internet_mediated_part <- 
  itanes_2013$d72_6 %in% c(1,2,3)
itanes_2013$internet_mediated_part[itanes_2013$d72_6 %in% c(5,6)] <- 
  NA

# TRUST 
itanes_2013$d15_1[itanes_2013$d15_1 > 4] <- NA
itanes_2013$d15_2[itanes_2013$d15_2 > 4] <- NA
itanes_2013$nat_pol_trust <- 
  rescale(as.numeric(itanes_2013$d15_1) + as.numeric(itanes_2013$d15_2),
          to = c(0,10), from = c(2,8))

weighted.mean(itanes_2013$nat_pol_trust, weight = weight, na.rm = T)

itanes_2018$nat_pol_trust <-
  rescale(as.numeric(itanes_2018$it_post_1_1_slice) + 
            as.numeric(itanes_2018$it_post_1_2_slice),
          to = c(0,10), from = c(2,22))

mean(itanes_2018$nat_pol_trust, na.rm  = T)

recodeIntPerTrust2013 <- function(x) {
  x <- as.numeric(x)
  x[x %in% c(8,9)] <- NA
  return(x)
}

itanes_2013$interpers_trust <-
  rescale(rowMeans(cbind(
    recodeIntPerTrust2013(itanes_2013$d18),
    recodeIntPerTrust2013(itanes_2013$d19),
    recodeIntPerTrust2013(itanes_2013$d20)),
    na.rm = T),
          to = c(0,10), from = c(1,7))

recodeIntPerTrust2018 <- function(x) {
  x <- as.numeric(x)
  x[x %in% c(12,13)] <- NA
  return(x)
}

itanes_2018$interpers_trust <-
    rescale(rowMeans(cbind(
    recodeIntPerTrust2018(itanes_2018$it_5_1_slice),
    recodeIntPerTrust2018(itanes_2018$it_6_1_slice)),
    na.rm = T),
          to = c(0,10), from = c(1,11))

# Association
itanes_2013$member <- 
  itanes_2013$d78_1 == 1 |
  itanes_2013$d78_5 == 1 |
  itanes_2013$d78_6 == 1 |
  itanes_2013$d78_7 == 1

# Political interest
itanes_2013$pol_interest <- 0
itanes_2013$pol_interest[itanes_2013$d27 == 3] <- 1
itanes_2013$pol_interest[itanes_2013$d27 == 2] <- 2
itanes_2013$pol_interest[itanes_2013$d27 == 1] <- 3
itanes_2013$pol_interest[itanes_2013$d27 == 6] <- NA  
itanes_2013$pol_interest <-
  rescale(itanes_2013$pol_interest, 
          to = c(0,10))

itanes_2018$pol_interest <- 0
itanes_2018$pol_interest[itanes_2018$d5_1 == 3] <- 1
itanes_2018$pol_interest[itanes_2018$d5_1 == 2] <- 2
itanes_2018$pol_interest[itanes_2018$d5_1 == 1] <- 3
itanes_2018$pol_interest[itanes_2018$d5_1 == 5] <- NA  
itanes_2018$pol_interest <-
  rescale(itanes_2018$pol_interest, 
          to = c(0,10))

# Use interne
itanes_2013$use_facebook <- 0
itanes_2013$use_facebook[itanes_2013$d73 == 3] <- 1
itanes_2013$use_facebook[itanes_2013$d73 == 2] <- 2
itanes_2013$use_facebook[itanes_2013$d73 == 1] <- 3
itanes_2013$use_facebook[itanes_2013$d73 == 6] <- NA

# Political talk
itanes_2013$poltalk_freq <- as.numeric(itanes_2013$d28)
itanes_2013$poltalk_freq[itanes_2013$d28 %in% c(6,7)] <- NA

# Demographics 
itanes_2013$female <- itanes_2013$sex == 2
itanes_2013$age <- itanes_2013$d1

itanes_2018$age <- as.numeric(itanes_2018$EtaCl)

# Ecological
itanes_2013$PRO_COM_T <-
  as.character(comuni_sp_data$PRO_COM_T)[match(itanes_2013$comune, 
                                 toupper(as.character(comuni_sp_data$COMUNE)))]
itanes_2013$PRO_COM_T[itanes_2013$comune == "SAINT VINCENT"] <- 
  "007065"
itanes_2013$PRO_COM_T[itanes_2013$comune == "ALBANO S ALESSANDRO"] <- 
  "016003"
itanes_2013$PRO_COM_T[itanes_2013$comune == "TOSCOLANO MADERNO"] <- 
  "017187"
itanes_2013$PRO_COM_T[itanes_2013$comune == "CAERANO SAN MARCO"] <- 
  "026006"
itanes_2013$PRO_COM_T[itanes_2013$comune == "MARANO DI VALPOLICEL"] <- 
  "023046"
itanes_2013$PRO_COM_T[itanes_2013$comune == "REGGIO EMILIA"] <- 
  "035033"
itanes_2013$PRO_COM_T[itanes_2013$comune == "PENNA S  ANDREA"] <- 
  "067033"
itanes_2013$PRO_COM_T[itanes_2013$comune == "CASTELLAMMARE DI STA"] <- 
  "063024"
itanes_2013$PRO_COM_T[itanes_2013$comune == "CAVA DEI TIRRENI"] <- 
  "065037"

itanes_2013 <- 
  merge(itanes_2013, 
        comuni_sp_data %>% 
          dplyr::select(PRO_COM_T, 
                        M5S_2013, M5S_2013_neigh_avg,
                        meetup_treated_2013,
                        meetup2014_19km_buff_freq_30days,
                        unemployment, partwf_perc, ANPIV_perc, 
                        pop_density, dist_10k_ppl,
                        ANPIV_perc, ref16_turnout, income2016_procapite), 
        by = 'PRO_COM_T')

itanes_2018$regint <-
  as.numeric(itanes_2018$regio1)

itanes_2018 <-
  merge(itanes_2018, 
        reg_long_df %>%
          filter(year == 2016) %>%
          dplyr::select(COD_REG, unemployment, employment,
                        n_meetup_events_per_10k,
                 ANPIV_perc, ref16_turnout, income),
        by.x = 'regint', by.y = 'COD_REG')

itanes_2013$employment <- itanes_2013$partwf_perc

itanes_2013$income_norm <- 
  rescale(itanes_2013$income2016_procapite, 
          to = c(0,10))
itanes_2018$income_norm <- 
  rescale(itanes_2018$income, 
          to = c(0,10))

itanes_2013_logit_1 <- 
  glm(M5S ~ 
        unemployed + 
        eco_stress + eco_perception +
        unemployment + employment +
        interpers_trust + 
        info_internet  + info_personal + info_tv + info_written +
        meetup_treated_2013 + nat_pol_trust +
        internet_mediated_part +
        M5S_2013 + M5S_2013_neigh_avg +
        ANPIV_perc + log(pop_density) + log(dist_10k_ppl) + ref16_turnout +
      income_norm + poltalk_freq,
      data = itanes_2013, family = quasibinomial, 
      weight = weight)

itanes_2013_logit_2 <- 
  glm(M5S ~ 
        unemployed + 
        eco_stress + eco_perception +
        unemployment + employment +
        interpers_trust + 
        info_internet  + info_personal + info_tv + info_written +
        meetup_treated_2013 + nat_pol_trust +
        internet_mediated_part +
        M5S_2013 + M5S_2013_neigh_avg +
        ANPIV_perc + log(pop_density) + log(dist_10k_ppl) + ref16_turnout +
      income_norm + poltalk_freq + age,
      data = itanes_2013, family = quasibinomial, 
      weight = weight)

itanes_2018_logit_1 <- 
  glm(M5S ~ 
        unemployed + 
        eco_perception +
        unemployment + employment +
        interpers_trust + 
        info_internet_news + info_social_media + 
        info_personal + info_tv + info_written +
        nat_pol_trust +
        ANPIV_perc +
        AmpQ + income_norm,
      data = itanes_2018, family = quasibinomial)

itanes_2018_logit_2 <- 
  glm(M5S ~ 
        unemployed + 
        eco_perception +
        unemployment + employment +
        interpers_trust + 
        info_internet_news + info_social_media + 
        info_personal + info_tv + info_written +
        nat_pol_trust +
        ANPIV_perc +
        AmpQ + income_norm + age,
      data = itanes_2018, family = quasibinomial)

itanes_2013$eco_stress_inv <- 
  rescale(itanes_2013$eco_stress, to = c(10,0))

itanes_2013_lm_2 <- 
  lm(poltalk_freq ~ 
       meetup_treated_2013 * nat_pol_trust +
       clsness_M5S  +
       eco_stress_inv +
       ANPIV_perc + ref16_turnout +
       pol_interest + degree + unemployed + 
       female + age
      + log(pop_density) + log(dist_10k_ppl)
     ,
     data = itanes_2013[itanes_2013$M5S == 1,], 
     weight = weight)

itanes_2013_lm_3 <- 
  lm(poltalk_freq ~ 
       meetup_treated_2013 * eco_stress_inv +
       clsness_M5S +
       nat_pol_trust +
       ANPIV_perc + ref16_turnout +
       pol_interest + degree + unemployed + 
       female + age
      + log(pop_density) + log(dist_10k_ppl)
     ,
     data = itanes_2013[itanes_2013$M5S == 1,], 
     weight = weight)

itanes_2013_lm_8 <- 
  lm(poltalk_freq ~ 
       meetup_treated_2013 + nat_pol_trust +
       clsness_M5S  +
       eco_stress_inv +
       ANPIV_perc + ref16_turnout +
       pol_interest + degree + unemployed + 
       female + age
      + log(pop_density) + log(dist_10k_ppl)
     ,
     data = itanes_2013[itanes_2013$M5S == 1,], 
     weight = weight)

itanes_2013_lm_4 <- 
  lm(poltalk_freq ~ 
       meetup_treated_2013 * nat_pol_trust * internet_mediated_part +
       clsness_M5S  +
       eco_stress +
       ANPIV_perc + ref16_turnout +
       pol_interest + degree + unemployed + 
       female + age
      + log(pop_density) + log(dist_10k_ppl)
     ,
     data = itanes_2013[itanes_2013$M5S == 1,], 
     weight = weight)

itanes_2013_lm_5 <- 
  lm(pol_part ~ 
       nat_pol_trust * clsness_M5S +
       eco_stress +
       # ANPIV_perc + ref16_turnout +
       member + interpers_trust +
       pol_interest + degree + unemployed + 
       female + age
     + log(pop_density) + log(dist_10k_ppl)
     ,
     data = itanes_2013, 
     weight = weight)

itanes_2013_lm_6 <- 
  lm(pol_part ~ 
       nat_pol_trust +
       eco_stress * clsness_M5S +
       # ANPIV_perc + ref16_turnout +
       member + interpers_trust +
       pol_interest + degree + unemployed + 
       female + age
     + log(pop_density) + log(dist_10k_ppl)
     ,
     data = itanes_2013, 
     weight = weight)

itanes_2013_lm_7 <- 
  lm(pol_part ~ 
       nat_pol_trust +
       eco_stress + clsness_M5S +
       # ANPIV_perc + ref16_turnout +
       member + interpers_trust +
       pol_interest + degree + unemployed + 
       female + age
     + log(pop_density) + log(dist_10k_ppl)
     ,
     data = itanes_2013, 
     weight = weight)

baseline_df <- 
  data.frame(nat_pol_trust = 0:10, 
             ANPIV_perc = 
               median(itanes_2013$ANPIV_perc[itanes_2013$M5S == 1]),
             ref16_turnout = 
               median(itanes_2013$ref16_turnout[itanes_2013$M5S == 1]),
             pol_interest = 
               median(itanes_2013$pol_interest[itanes_2013$M5S == 1]),
             degree = 
               median(itanes_2013$degree[itanes_2013$M5S == 1]),
             unemployed = 
               median(itanes_2013$unemployed[itanes_2013$M5S == 1]),
             female = 
               TRUE,
             age = 
               median(itanes_2013$age[itanes_2013$M5S == 1]),
             pop_density = 
               median(itanes_2013$pop_density[itanes_2013$M5S == 1]),
             dist_10k_ppl = 
               median(itanes_2013$dist_10k_ppl[itanes_2013$M5S == 1]),
             clsness_M5S = 
               median(itanes_2013$clsness_M5S[itanes_2013$M5S == 1],
                      na.rm = T),
             eco_stress_inv = 
               median(itanes_2013$eco_stress_inv[itanes_2013$M5S == 1],
                      na.rm = T))

itanes_2013_lm_2_predTRUE <-
  predict(itanes_2013_lm_2, 
          newdata = cbind(baseline_df, meetup_treated_2013=TRUE),
          se = TRUE)

itanes_2013_lm_2_predFALSE <-
  predict(itanes_2013_lm_2, 
          newdata = cbind(baseline_df, meetup_treated_2013=FALSE),
          se = TRUE)

itanes_2013_lm_2_pred.df <-
  rbind(with(itanes_2013_lm_2_predTRUE,
             data.frame(nat_pol_trust = 0:10,
                        meetup_treated_2013 = TRUE,     
                        fit,
                        upper_err = fit + se.fit,
                        lower_err = fit - se.fit)),
        with(itanes_2013_lm_2_predFALSE,
             data.frame(nat_pol_trust = 0:10,
                        meetup_treated_2013 = FALSE,       
                        fit,
                        upper_err = fit + se.fit,
                        lower_err = fit - se.fit))
  )

predict_itanes_logit_2013 <- 
  itanes_2013 %>%
  dplyr::select(unemployed, 
                eco_stress, eco_perception,
                unemployment, employment,
                interpers_trust, 
                info_internet , info_personal, info_tv, info_written,
                meetup_treated_2013, nat_pol_trust,
                internet_mediated_part,
                M5S_2013, M5S_2013_neigh_avg,
                ANPIV_perc, pop_density, dist_10k_ppl, ref16_turnout,
                income_norm, poltalk_freq) %>%
  summarize_all(median, na.rm = T)
predict_itanes_logit_2013$unemployed <- 
  as.logical(predict_itanes_logit_2013$unemployed)
predict_itanes_logit_2013$info_internet <- 
  as.logical(predict_itanes_logit_2013$info_internet)
predict_itanes_logit_2013$info_personal <- 
  as.logical(predict_itanes_logit_2013$info_personal)
predict_itanes_logit_2013$info_internet <- 
  as.logical(predict_itanes_logit_2013$info_internet)
predict_itanes_logit_2013$info_tv <- 
  as.logical(predict_itanes_logit_2013$info_tv)
predict_itanes_logit_2013$info_written <- 
  as.logical(predict_itanes_logit_2013$info_written)
predict_itanes_logit_2013$meetup_treated_2013 <- 
  as.logical(predict_itanes_logit_2013$meetup_treated_2013)

prop.table(table(itanes_2013$M5S))

predict(itanes_2013_logit_1, 
        predict_itanes_logit_2013[rep(1, 4),] %>%
          mutate(nat_pol_trust = c(0, 2.705, 5, 10)), 
        type = "response")

predict(itanes_2013_logit_1, 
        predict_itanes_logit_2013[rep(1, 2),] %>%
          mutate(internet_mediated_part = c(TRUE, FALSE)), 
        type = "response")

predict(itanes_2013_logit_1, 
        predict_itanes_logit_2013[rep(1, 2),] %>%
          mutate(poltalk_freq = c(1, 5)), 
        type = "response")
```

### Table 4

```{r}
  
do.call(
  stargazer, 
  c(list(itanes_2013_logit_1, 
         itanes_2018_logit_1,
         # type = 'latex', out = "table/itanes_2013-18_logit.tex",
         type = 'html', out = "table/itanes_2013-18_logit.html",
         dep.var.labels = "M5S vote (binary)",
         column.labels = c("2013 (B)", "2013 (B)", "2018 (B)", "2018 (B)"),
         omit = c('unemployment', 'employment', 
                  "interpers_trust", "info_personal",
                  "info_tv", "info_written", "M5S_2013", 
                  "M5S_2013_neigh_avg", "pop_density", "dist_10k_ppl",
                  "ref16_turnout", 'Constant', "ANPIV_perc", "income_norm", "AmpQ",
                  "info_internet_news"),
         order = c("nat_pol_trust", 
                   "eco_stress", "eco_perception","unemployed",
                   "info_internet", "info_social_media",
                   "internet_mediated_part",
                   "poltalk_freq",
                   "meetup_treated_2013"),
         covariate.labels = var_dict_labels[
           c("nat_pol_trust",
             "eco_stress", "eco_perception","unemployed",
             "info_internet", "info_social_media",
             "internet_mediated_part",
             "poltalk_freq",
             "meetup_treated_2013")]),
    sg_param_intext))
```


```{r}
# Issue salience

## 2013
itanes_salience_firstone_2013 <-
  data.frame(M5S  = itanes_2013$M5S == 1, 
             `Employment policies` = 
               as_factor(itanes_2013$d11_rec) == "Employment policies",
             `Economic development` = 
               as_factor(itanes_2013$d11_rec) == "Economic development",
             `Political ethics` = 
               as_factor(itanes_2013$d11_rec) == "Political ethics",
             `Taxes` = 
               as_factor(itanes_2013$d11_rec) == "Taxes",
             `Economic insecurity` = 
               as_factor(itanes_2013$d11_rec) == "Economic insecurity",
             `Business aids` = 
               as_factor(itanes_2013$d11_rec) == "Business aids",
             `Public administration/Institutional reforms` = 
               as_factor(itanes_2013$d11_rec) == "Public administration/Institutional reforms" | 
               as_factor(itanes_2013$d12_rec) == "Public administration/Institutional reforms",
             `Legality` = 
               as_factor(itanes_2013$d11_rec) == "Legality"
  )

itanes_salience_firsttwo_2013 <-
  data.frame(M5S  = itanes_2013$M5S == 1, 
             `Employment policies` = 
               as_factor(itanes_2013$d11_rec) == "Employment policies" | 
               as_factor(itanes_2013$d12_rec) == "Employment policies",
             `Economic development` = 
               as_factor(itanes_2013$d11_rec) == "Economic development" | 
               as_factor(itanes_2013$d12_rec) == "Economic development",
             `Political ethics` = 
               as_factor(itanes_2013$d11_rec) == "Political ethics" | 
               as_factor(itanes_2013$d12_rec) == "Political ethics",
             `Taxes` = 
               as_factor(itanes_2013$d11_rec) == "Taxes" | 
               as_factor(itanes_2013$d12_rec) == "Taxes",
             `Economic insecurity` = 
               as_factor(itanes_2013$d11_rec) == "Economic insecurity" | 
               as_factor(itanes_2013$d12_rec) == "Economic insecurity",
             `Business aids` = 
               as_factor(itanes_2013$d11_rec) == "Business aids" | 
               as_factor(itanes_2013$d12_rec) == "Business aids",
             `Public administration/Institutional reforms` = 
               as_factor(itanes_2013$d11_rec) == "Public administration/Institutional reforms" | 
               as_factor(itanes_2013$d12_rec) == "Public administration/Institutional reforms",
             `Legality` = 
               as_factor(itanes_2013$d11_rec) == "Legality" | 
               as_factor(itanes_2013$d12_rec) == "Legality"
  )
  

itanes_salience_firstone_stats_2013 <- data.frame()

for (this_colname in names(itanes_salience_firstone_2013)[2:9]) {
  this_df <-
    itanes_salience_firstone_2013 %>%
    dplyr::group_by(M5S) %>%
    dplyr::summarize(mean = mean(!!sym(this_colname), na.rm = T),
                     sd = sd(!!sym(this_colname), na.rm = T),
                     se = sd(!!sym(this_colname), na.rm = T) / 
                       sqrt(sum(!is.na(!!sym(this_colname)))))
  this_df$issue <- this_colname
  itanes_salience_firstone_stats_2013 <-
    rbind(itanes_salience_firstone_stats_2013, 
          this_df)
}

itanes_salience_firstone_stats_2013$issue <-
  gsub("\\.", " ", itanes_salience_firstone_stats_2013$issue)
itanes_salience_firstone_stats_2013$issue[grepl("Institutional reforms",
                                                itanes_salience_firstone_stats_2013$issue)] <-
  "Govt./Institutional reforms"

itanes_salience_firsttwo_stats_2013 <- data.frame()

for (this_colname in names(itanes_salience_firsttwo_2013)[2:9]) {
  this_df <-
    itanes_salience_firsttwo_2013 %>%
    dplyr::group_by(M5S) %>%
    dplyr::summarize(mean = mean(!!sym(this_colname), na.rm = T),
                     sd = sd(!!sym(this_colname), na.rm = T),
                     se = sd(!!sym(this_colname), na.rm = T) / 
                       sqrt(sum(!is.na(!!sym(this_colname)))))
  this_df$issue <- this_colname
  itanes_salience_firsttwo_stats_2013 <-
    rbind(itanes_salience_firsttwo_stats_2013, 
          this_df)
}

itanes_salience_firsttwo_stats_2013$issue <-
  gsub("\\.", " ", itanes_salience_firsttwo_stats_2013$issue)
itanes_salience_firsttwo_stats_2013$issue[grepl("Institutional reforms",
                                                itanes_salience_firsttwo_stats_2013$issue)] <-
  "Govt./Institutional reforms"


## 2018

itanes_salience_firstone_2018 <-
  data.frame(M5S  = itanes_2018$M5S == 1, 
             `Taxes` = 
               as_factor(itanes_2018$it_3) == "Le tasse",
                          `Unemployment` = 
               as_factor(itanes_2018$it_3) == "La disoccupazione",
                                       `Public debt` = 
               as_factor(itanes_2018$it_3) == "Il debito pubblico",
                                       `EU economic policies` = 
               as_factor(itanes_2018$it_3) == "Le politiche economiche dell’Europa",
                                       `Tax evasion` = 
               as_factor(itanes_2018$it_3) == "L’evasione fiscale",
                                       `Political corruption` = 
               as_factor(itanes_2018$it_3) == "La corruzione politica",
                                       `Immigration` = 
               as_factor(itanes_2018$it_3) == "L’immigrazione",
                                       `Crime` = 
               as_factor(itanes_2018$it_3) == "La criminalità",
                                       `Economic stagnation` = 
               as_factor(itanes_2018$it_3) == "La mancanza di crescita",
                                       `Environment protection` = 
               as_factor(itanes_2018$it_3) == "La difesa dell’ambiente"
  )


itanes_salience_firsttwo_2018 <-
  data.frame(M5S  = itanes_2018$M5S == 1, 
             `Taxes` = 
               as_factor(itanes_2018$it_3) == "Le tasse" | 
               as_factor(itanes_2018$it_3_1) == "Le tasse",
                          `Unemployment` = 
               as_factor(itanes_2018$it_3) == "La disoccupazione" | 
               as_factor(itanes_2018$it_3_1) == "La disoccupazione",
                                       `Public debt` = 
               as_factor(itanes_2018$it_3) == "Il debito pubblico" | 
               as_factor(itanes_2018$it_3_1) == "Il debito pubblico",
                                       `EU economic policies` = 
               as_factor(itanes_2018$it_3) == "Le politiche economiche dell’Europa" | 
               as_factor(itanes_2018$it_3_1) == "Le politiche economiche dell’Europa",
                                       `Tax evasion` = 
               as_factor(itanes_2018$it_3) == "L’evasione fiscale" | 
               as_factor(itanes_2018$it_3_1) == "L’evasione fiscale",
                                       `Political corruption` = 
               as_factor(itanes_2018$it_3) == "La corruzione politica" | 
               as_factor(itanes_2018$it_3_1) == "La corruzione politica",
                                       `Immigration` = 
               as_factor(itanes_2018$it_3) == "L’immigrazione" | 
               as_factor(itanes_2018$it_3_1) == "L’immigrazione",
                                       `Crime` = 
               as_factor(itanes_2018$it_3) == "La criminalità" | 
               as_factor(itanes_2018$it_3_1) == "La criminalità",
                                       `Economic stagnation` = 
               as_factor(itanes_2018$it_3) == "La mancanza di crescita" | 
               as_factor(itanes_2018$it_3_1) == "La mancanza di crescita",
                                       `Environment protection` = 
               as_factor(itanes_2018$it_3) == "La difesa dell’ambiente" | 
               as_factor(itanes_2018$it_3_1) == "La difesa dell’ambiente"
  )


itanes_salience_firstone_stats_2018 <- data.frame()

for (this_colname in names(itanes_salience_firstone_2018)[2:9]) {
  this_df <-
    itanes_salience_firstone_2018 %>%
    dplyr::group_by(M5S) %>%
    dplyr::summarize(mean = mean(!!sym(this_colname), na.rm = T),
                     sd = sd(!!sym(this_colname), na.rm = T),
                     se = sd(!!sym(this_colname), na.rm = T) / 
                       sqrt(sum(!is.na(!!sym(this_colname)))))
  this_df$issue <- this_colname
  itanes_salience_firstone_stats_2018 <-
    rbind(itanes_salience_firstone_stats_2018, 
          this_df)
}

itanes_salience_firsttwo_stats_2018 <- data.frame()

for (this_colname in names(itanes_salience_firsttwo_2018)[2:9]) {
  this_df <-
    itanes_salience_firsttwo_2018 %>%
    dplyr::group_by(M5S) %>%
    dplyr::summarize(mean = mean(!!sym(this_colname), na.rm = T),
                     sd = sd(!!sym(this_colname), na.rm = T),
                     se = sd(!!sym(this_colname), na.rm = T) / 
                       sqrt(sum(!is.na(!!sym(this_colname)))))
  this_df$issue <- this_colname
  itanes_salience_firsttwo_stats_2018 <-
    rbind(itanes_salience_firsttwo_stats_2018, 
          this_df)
}


itanes_salience_final_stats_2013_2018 <-
  rbind(itanes_salience_firstone_stats_2013 %>%
          mutate(year = '2013', type = "As 1st issue"), 
        itanes_salience_firstone_stats_2018 %>%
          mutate(year = '2018', type = "As 1st issue"),
        itanes_salience_firsttwo_stats_2013 %>%
          mutate(year = '2013', type = "As 1st or 2nd issue"),
        itanes_salience_firsttwo_stats_2018 %>%
          mutate(year = '2018', type = "As 1st or 2nd issue"))

itanes_salience_final_stats_2013_2018$issue <-
  gsub("\\.", " ", itanes_salience_final_stats_2013_2018$issue)

itanes_salience_final_stats_2013_2018 <-
  itanes_salience_final_stats_2013_2018 %>%
  dplyr::filter(!is.na(M5S)) %>%
  dplyr::group_by(issue, year) %>%
  dplyr::mutate(diff = mean[M5S == TRUE & 
                              type == "As 1st issue"] - 
                  mean[M5S == FALSE &
                         type == "As 1st issue"],
                se_1st = se[M5S == TRUE &
                              type == "As 1st issue"])

itanes_salience_final_stats_2013_2018$issue <-
  paste0(itanes_salience_final_stats_2013_2018$issue,
         "\nDiff.=", 
         round(itanes_salience_final_stats_2013_2018$diff*100, 2), 
         "pp, SE=",
         round(itanes_salience_final_stats_2013_2018$se_1st*100, 2),"pp")

itanes_salience_final_stats_2013_2018$issue <-
  factor(itanes_salience_final_stats_2013_2018$issue,
         levels = unique(itanes_salience_final_stats_2013_2018$issue[
           order(itanes_salience_final_stats_2013_2018$diff, 
                 decreasing = TRUE)
         ]))

itanes_salience_final_stats_2013_2018$M5S <-
  factor(itanes_salience_final_stats_2013_2018$M5S,
         levels = c(TRUE, FALSE))
```

```{r itanes_salience_2013, fig.width =  7, fig.height = 7}
ggplot(itanes_salience_final_stats_2013_2018 %>%
         filter(year == '2013'), 
       aes(x=type, y=mean, fill = M5S,
           ymin=mean-se, ymax=mean+se)) +
  geom_bar(stat='identity', position = 'dodge2') +
  geom_errorbar(
              width=.2, position=position_dodge(.9)) +
  facet_wrap(.~issue, scales = "free_y") +
  scale_y_continuous(label = percent) +
  scale_fill_manual(labels = c("Yes", "No"), values = two_colours) +
  theme_bw() +
  labs(x = NULL, y = NULL, fill = "Voted M5S in 2013") +
  theme(legend.position = c(.85,.15))
```

```{r itanes_salience_2018, fig.width =  7, fig.height = 7}
ggplot(itanes_salience_final_stats_2013_2018 %>%
         filter(!is.na(M5S) & 
                  year == '2018'), 
       aes(x=type, y=mean, fill = M5S,
           ymin=mean-se, ymax=mean+se)) +
  geom_bar(stat='identity', position = 'dodge2') +
  geom_errorbar(
              width=.2, position=position_dodge(.9)) +
  facet_wrap(.~issue, scales = "free_y") +
  scale_y_continuous(label = percent) +
  scale_fill_manual(labels = c("Yes", "No"), values = two_colours) +
  theme_bw() +
  labs(x = NULL, y = NULL, fill = "Voted M5S in 2018") +
  theme(legend.position = c(.85,.15))
```

### Figure 5

```{r itanes_interaction, fig.width =  6, fig.height = 8}
library(egg)
ggarrange(
  ggplot(itanes_2013_lm_2_pred.df,
         aes(x = nat_pol_trust, fill = meetup_treated_2013,
             y = fit,
             ymax = upper_err,
             ymin = lower_err)) +
    geom_line() +
    geom_ribbon(alpha = .5) +
    geom_vline(xintercept = quantile(itanes_2013$nat_pol_trust[itanes_2013$M5S ==
                                                               1], probs = .5, na.rm = T),
               colour = 'red', linetype=3) +
    labs(x = var_dict_labels[['nat_pol_trust']],
         y = var_dict_labels[['poltalk_freq']],
         fill = var_dict_labels[['meetup_treated_2013']]) +
    scale_fill_manual(values = two_colours) +
    theme_bw() +
    theme(legend.position = c(.45,.23), 
          legend.background = element_rect(fill = "transparent")),
  ggplot(itanes_2013[itanes_2013$M5S == 1,], 
         aes(y=(..count..)/sum(..count..), x = nat_pol_trust)) +
    geom_bar() + theme_bw() +
    scale_y_continuous(labels = percent_format(accuracy = 1)) +
    labs(x = var_dict_labels[['nat_pol_trust']], y = NULL),
  ncol = 1)
```

### Table 5

```{r, results='asis'}
do.call(stargazer,
        c(list(itanes_2013_lm_2, itanes_2013_lm_3,
               # type = 'latex', out = "table/itanes_2013_lm.tex",
               type = 'html', out = "table/itanes_2013_lm.html",
               dep.var.labels = var_dict_labels['poltalk_freq'],
               column.labels = c("1 (B)", "2 (B)"),
               omit = c('clsness_M5S', 'pop_density',
                        'dist_10k_ppl', 'ref16_turnout', 'age'),
               order = c("^meetup_treated_2013TRUE$",
                         "^nat_pol_trust$",
                         "^eco_stress_inv$", 
                         "pol_interest", 
                         "degree", "unemployed", "female",
                         "ANPIV_perc",
                         "meetup_treated_2013TRUE:nat_pol_trust",
                         "meetup_treated_2013TRUE:eco_stress_inv"),
               covariate.labels = var_dict_labels[
                 c("meetup_treated_2013",
                   "nat_pol_trust",
                   "eco_stress_inv",
                   "pol_interest",
                   "degree", "unemployed", "female",
                   "ANPIV_perc",
                   "meetup_treated_2013TRUE:nat_pol_trust",
                   "meetup_treated_2013TRUE:eco_stress_inv")]),
          sg_param_intext))
```

### Table 6 (Supplemental online material)

```{r, results='asis'}
do.call(
  stargazer,
  c(list(itanes_2013_lm_2, itanes_2013_lm_3, itanes_2013_lm_8,
         # type = 'latex', out = "table/itanes_2013_lm_full.tex",
         type = 'html', out = "table/itanes_2013_lm_full.html",
         dep.var.labels = var_dict_labels['poltalk_freq'],
         column.labels = c("(B)", "(B)", "(B)"),
         order = c("^meetup_treated_2013TRUE$",
                   "^nat_pol_trust$",
                   "^eco_stress_inv$", 
                   "pol_interest", 
                   "degree", "unemployed", "female",
                   "ANPIV_perc",
                   "meetup_treated_2013TRUE:nat_pol_trust",
                   "meetup_treated_2013TRUE:eco_stress_inv")
         ,
         covariate.labels = var_dict_labels[
           c("meetup_treated_2013",
             "nat_pol_trust",
             "eco_stress_inv",
             "pol_interest",
             "degree", "unemployed", "female",
             "ANPIV_perc",
             "meetup_treated_2013TRUE:nat_pol_trust",
             "meetup_treated_2013TRUE:eco_stress_inv",
             "clsness_M5S", "ref16_turnout", "age", "log(pop_density)",
             "log(dist_10k_ppl)")]
         ),
    sg_param_appendix))
```

### Table 2

```{r}
itanes_2013_lm_5_stdcoef <- 
  lm.beta(itanes_2013_lm_5)
do.call(
  stargazer,
  c(list(itanes_2013_lm_5,
          # type = 'latex', out = "table/itanes_2013_lm_pol_part.tex",
          type = 'html', out = "table/itanes_2013_lm_pol_part.html",
         dep.var.labels = var_dict_labels['pol_part'],
          column.labels = "(Beta)",
          omit = c('pop_density', 'female',
                   'dist_10k_ppl', 'age', 'nat_pol_trust'),
          coef = list(itanes_2013_lm_5_stdcoef$standardized.coefficients),
          p = list(summary(itanes_2013_lm_5)$coefficients[,4]),
          covariate.labels = c("Closeness to M5S | Trust in nat. inst. = 0",
                               var_dict_labels[
            c("eco_stress",
              "member", "interpers_trust",
              "pol_interest",
                    "degree", "unemployed",
                    "nat_pol_trust:clsness_M5S")])),
    sg_param_intext))
```

### Table 8 (Supplemental online material)

```{r}
itanes_2013_lm_5_stdcoef <- 
  lm.beta(itanes_2013_lm_5)
itanes_2013_lm_6_stdcoef <- 
  lm.beta(itanes_2013_lm_6)
itanes_2013_lm_7_stdcoef <- 
  lm.beta(itanes_2013_lm_7)

do.call(
  stargazer,
  c(list(itanes_2013_lm_5, itanes_2013_lm_5,
         itanes_2013_lm_6, itanes_2013_lm_6,
         itanes_2013_lm_7, itanes_2013_lm_7,
         coef = list(itanes_2013_lm_5$coefficients,
                     itanes_2013_lm_5_stdcoef$standardized.coefficients,
                     itanes_2013_lm_6$coefficients,
                     itanes_2013_lm_6_stdcoef$standardized.coefficients,
                     itanes_2013_lm_7$coefficients,
                     itanes_2013_lm_7_stdcoef$standardized.coefficients),
         p = list(summary(itanes_2013_lm_5)$coefficients[,4],
                  summary(itanes_2013_lm_5)$coefficients[,4],
                  summary(itanes_2013_lm_6)$coefficients[,4],
                  summary(itanes_2013_lm_6)$coefficients[,4],
                  summary(itanes_2013_lm_7)$coefficients[,4],
                  summary(itanes_2013_lm_7)$coefficients[,4]),
         column.labels = rep(c("(B)", "(Beta)"), 3),
         type = 'latex', out = "table/itanes_2013_lm_pol_part_full.tex",
         dep.var.labels = var_dict_labels['pol_part'],
         covariate.labels = var_dict_labels[
           c("nat_pol_trust", "clsness_M5S", "eco_stress",
             "member", "interpers_trust",
             "pol_interest",
             "degree", "unemployed", "female", "age",
             "log(pop_density)", "log(dist_10k_ppl)",
             "nat_pol_trust:clsness_M5S", "eco_stress:clsness_M5S")]),
    sg_param_appendix))
```

### Table 3 (Supplemental online material)

```{r, results = 'asis'}
do.call(
stargazer, c(list(itanes_2013_logit_1, itanes_2013_logit_2, 
                  itanes_2018_logit_1, itanes_2018_logit_2,
          type = 'latex', out = "table/itanes_2013-18_logit_full.tex",
          # type = 'html', out = "table/itanes_2013-18_logit_full.html",
          dep.var.labels = "M5S vote",
          column.labels = c("2013 (bin.)", "2013 (bin.)", 
                            "2018 (bin.)", "2018 (bin.)")
          ,
          covariate.labels = var_dict_labels[
            c('unemployed','eco_stress', 'eco_perception',
              'unemployment', 'employment', 'interpers_trust',
              'info_internet', 'info_internet_news', 'info_social_media',
              'info_personal', 'info_tv', 'info_written', 'meetup_treated_2013',
              'nat_pol_trust', 'internet_mediated_part', 'M5S_2013',
              'M5S_2013_neigh_avg',
              'ANPIV_perc', 'log(pop_density)', 'log(dist_10k_ppl)', 'ref16_turnout',
              'AmpQ', 'income_norm', 'poltalk_freq', 'age')]
          ),
          sg_param_appendix))
```

```{r}
stargazer(itanes_2013_lm_2, itanes_2013_lm_3, itanes_2013_lm_4,
          type = 'latex', out = "table/itanes_2013_lm_full.tex",
          dep.var.labels = var_dict_labels['poltalk_freq'],
          omit.stat = c('adj.rsq', 'f', 'ser'),
          order = c("^meetup_treated_2013TRUE$",
                    "^nat_pol_trust$",
                    "^eco_stress$",
                    "^internet_mediated_partTRUE$",
                    "^pol_interest$",
                    "clsness_M5S",
                    "degree", "unemployed", "female", "age",
                    "ANPIV_perc", "ref16_turnout",
                    "pop_density", "dist_10k_ppl",
                    "meetup_treated_2013TRUE:nat_pol_trust^",
                    "meetup_treated_2013TRUE:eco_stress^",
                    "meetup_treated_2013TRUE:internet_mediated_part^",
                    "nat_pol_trust:internet_mediated_part^",
                    "meetup_treated_2013TRUE:nat_pol_trust:internet_mediated_part"),
          covariate.labels = var_dict_labels[
            c("meetup_treated_2013",
              "nat_pol_trust",
              "eco_stress",
              "internet_mediated_part",
              "pol_interest",
              "clsness_M5S",
              "degree", "unemployed", "female", "age",
              "ANPIV_perc", "ref16_turnout",
              "log(pop_density)", "log(dist_10k_ppl)",
              "meetup_treated_2013TRUE:nat_pol_trust",
              "meetup_treated_2013TRUE:eco_stress",
              "meetup_treated_2013TRUE:internet_mediated_part",
              "nat_pol_trust:internet_mediated_part",
              "meetup_treated_2013TRUE:nat_pol_trust:internet_mediated_part")],
          report = "vc*",
          model.numbers = FALSE,
          model.names = FALSE,
          float = FALSE,
          style = 'ajps')
```

```{r}
library(haven)
dat_2014 <- 
  read_sav("secondary_data/misc_european_elections/ZA5160_v4-0-0.sav")
dat_2009 <- 
  read_sav("secondary_data/misc_european_elections/ZA5055_v1-1-1.sav")

parties_2014 <- unique(dat_2014$qp2_ees)
partie_2009 <- unique(dat_2009$q25)
new_parties <- parties_2014[!parties_2014 %in% partie_2009]
old_parties <- parties_2014[parties_2014 %in% partie_2009]
  
dat_2014$internet_daily <- dat_2014$d62_1 == 1

require(dplyr)
dat_2014_new_parties <- dat_2014 %>% filter(qp2_ees %in% new_parties)
dat_2014_old_parties <- dat_2014 %>% filter(qp2_ees %in% old_parties)

weighted.mean(dat_2014_new_parties$internet_daily, 
              w = dat_2014_new_parties$wex)

weighted.mean(dat_2014_old_parties$internet_daily, 
              w = dat_2014_old_parties$wex)

# Internet

parties <- 
  c(`Podemos` = 1724230, 
    `Citizens` = 1724410,
    `Five Star Movement` = 1380956, 
    `Democratic Left` = 1300225,
    `Independent Greeks` = 1300611,
    `Golden Dawn` = 1300710,
    `Alternative for Germany` = 1276621,
    `The British National Party` = 1826720)

dat_2014$new_party <- dat_2014$qp2_ees %in% parties

require(weights)

is_spain <- dat_2014$countrycode == 1724
is_greece <- dat_2014$countrycode == 1300
is_italy <- dat_2014$countrycode == 1380
is_germany <- dat_2014$countrycode == 1276
is_uk <- dat_2014$countrycode == 1826

res <- data.frame()

test_res <- 
  wtd.t.test(dat_2014$internet_daily[is_spain & dat_2014$new_party],
             dat_2014$internet_daily[is_spain & !dat_2014$new_party],
             weight = dat_2014$wex[is_spain & dat_2014$new_party],
             weighty = dat_2014$wex[is_spain & !dat_2014$new_party])

res <- rbind(res, 
             data.frame(country = "Spain",
                        party = "All new parties",
                        `vote %` = "34.6 (2015)",
                        `Internet daily %` = test_res$additional[2],
                        `Diff.` = test_res$additional[1],
                        `p` = test_res$coefficients[3],
                        `test` = test_res$test,
                        `source` = "EES 2014"))

test_res <- 
  wtd.t.test(dat_2014$internet_daily[is_spain & dat_2014$qp2_ees == parties['Podemos']],
             dat_2014$internet_daily[is_spain & !dat_2014$new_party],
             weight = dat_2014$wex[is_spain & dat_2014$qp2_ees == parties['Podemos']],
             weighty = dat_2014$wex[is_spain & !dat_2014$new_party])

res <- rbind(res, 
             data.frame(country = "Spain",
                        party = "Podemos",
                        `vote %` = "20.7 (2015)",
                        `Internet daily %` = test_res$additional[2],
                        `Diff.` = test_res$additional[1],
                        `p` = test_res$coefficients[3],
                        `test` = test_res$test,
                        `source` = "EES 2014"))

test_res <- 
  wtd.t.test(dat_2014$internet_daily[is_spain & dat_2014$qp2_ees == parties['Citizens']],
             dat_2014$internet_daily[is_spain & !dat_2014$new_party],
             weight = dat_2014$wex[is_spain & dat_2014$qp2_ees == parties['Citizens']],
             weighty = dat_2014$wex[is_spain & !dat_2014$new_party])

res <- rbind(res, 
             data.frame(country = "Spain",
                        party = "Ciudadanos",
                        `vote %` = "13.9 (2015)",
                        `Internet daily %` = test_res$additional[2],
                        `Diff.` = test_res$additional[1],
                        `p` = test_res$coefficients[3],
                        `test` = test_res$test,
                        `source` = "EES 2014"))

test_res <- 
  wtd.t.test(dat_2014$internet_daily[is_greece & dat_2014$qp2_ees %in% new_parties],
             dat_2014$internet_daily[is_greece & !dat_2014$new_party],
             weight = dat_2014$wex[is_greece & dat_2014$qp2_ees %in% new_parties],
             weighty = dat_2014$wex[is_greece & !dat_2014$new_party])

res <- rbind(res, 
             data.frame(country = "Greece",
                        party = "All new parties",
                        `vote %` = "25.4 (May 2012)",
                        `Internet daily %` = test_res$additional[2],
                        `Diff.` = test_res$additional[1],
                        `p` = test_res$coefficients[3],
                        `test` = test_res$test,
                        `source` = "EES 2014"))

test_res <- 
  wtd.t.test(dat_2014$internet_daily[is_greece & dat_2014$qp2_ees == parties['Independent Greeks']],
             dat_2014$internet_daily[is_greece & !dat_2014$new_party],
             weight = dat_2014$wex[is_greece & dat_2014$qp2_ees == parties['Independent Greeks']],
             weighty = dat_2014$wex[is_greece & !dat_2014$new_party])

res <- rbind(res, 
             data.frame(country = "Greece",
                        party = "Independent Greeks",
                        `vote %` = "10.6 (May 2012)",
                        `Internet daily %` = test_res$additional[2],
                        `Diff.` = test_res$additional[1],
                        `p` = test_res$coefficients[3],
                        `test` = test_res$test,
                        `source` = "EES 2014"))

test_res <- 
  wtd.t.test(dat_2014$internet_daily[is_greece & dat_2014$qp2_ees == parties['Democratic Left']],
             dat_2014$internet_daily[is_greece & !dat_2014$new_party],
             weight = dat_2014$wex[is_greece & dat_2014$qp2_ees == parties['Democratic Left']],
             weighty = dat_2014$wex[is_greece & !dat_2014$new_party])

res <- rbind(res, 
             data.frame(country = "Greece",
                        party = "Democratic Left",
                        `vote %` = "6.1 (May 2012)",
                        `Internet daily %` = test_res$additional[2],
                        `Diff.` = test_res$additional[1],
                        `p` = test_res$coefficients[3],
                        `test` = test_res$test,
                        `source` = "EES 2014"))

test_res <- 
  wtd.t.test(dat_2014$internet_daily[is_greece & dat_2014$qp2_ees == parties['Golden Dawn']],
             dat_2014$internet_daily[is_greece & !dat_2014$new_party],
             weight = dat_2014$wex[is_greece & dat_2014$qp2_ees == parties['Golden Dawn']],
             weighty = dat_2014$wex[is_greece & !dat_2014$new_party])

res <- rbind(res, 
             data.frame(country = "Greece",
                        party = "Golden Dawn",
                        `vote %` = "7 (May 2012)",
                        `Internet daily %` = test_res$additional[2],
                        `Diff.` = test_res$additional[1],
                        `p` = test_res$coefficients[3],
                        `test` = test_res$test,
                        `source` = "EES 2014"))

test_res <- 
  wtd.t.test(dat_2014$internet_daily[is_italy & dat_2014$new_party],
             dat_2014$internet_daily[is_italy & !dat_2014$new_party],
             weight = dat_2014$wex[is_italy & dat_2014$new_party],
             weighty = dat_2014$wex[is_italy & !dat_2014$new_party])

res <- rbind(res, 
             data.frame(country = "Italy",
                        party = "Five Star Movement",
                        `vote %` = "25.6 (2013)",
                        `Internet daily %` = test_res$additional[2],
                        `Diff.` = test_res$additional[1],
                        `p` = test_res$coefficients[3],
                        `test` = test_res$test,
                        `source` = "EES 2014"))

test_res <- 
  wtd.t.test(dat_2014$internet_daily[is_germany & dat_2014$qp2_ees == parties['Alternative for Germany']],
             dat_2014$internet_daily[is_germany & !dat_2014$new_party],
             weight = dat_2014$wex[is_germany & dat_2014$qp2_ees == parties['Alternative for Germany']],
             weighty = dat_2014$wex[is_germany & !dat_2014$new_party])

res <- rbind(res, 
             data.frame(country = "Germany",
                        party = "Alternative for Germany",
                        `vote %` = "4.7 (2013)",
                        `Internet daily %` = test_res$additional[2],
                        `Diff.` = test_res$additional[1],
                        `p` = test_res$coefficients[3],
                        `test` = test_res$test,
                        `source` = "EES 2014"))

# Iceland

dat_2013 <- 
  read_sav("secondary_data/misc_european_elections/icenes_2013_open_access_english_1release.sav")
dat_2013$internet_daily <- dat_2013$internet == 1

iceland_parties <- 
  c(`Bright Future` = 5,
    `Pirate Party` = 6, `Household Party` = 14,
    `Iceland Democratic Party` = 8, 
    `Right-Green People's Party` = 10, 
    `Rainbow` = 13)

dat_2013$new_party <- dat_2013$prtvote13 %in% iceland_parties

test_res <-
  t.test(dat_2013$internet_daily[dat_2013$new_party],
         dat_2013$internet_daily[!dat_2013$new_party])

res <- rbind(res, 
             data.frame(country = "Iceland",
                        party = "All new parties",
                        `vote %` = "21.7 (2013)",
                        `Internet daily %` = test_res$estimate[1],
                        `Diff.` = test_res$estimate[1] - test_res$estimate[2],
                        `p` = test_res$p.value,
                        `test` = "Student's t-Test",
                        `source` = "ICENES 2013"))

test_res <-
  t.test(dat_2013$internet_daily[dat_2013$prtvote13 == iceland_parties['Bright Future']],
       dat_2013$internet_daily[!dat_2013$new_party])

res <- rbind(res, 
             data.frame(country = "Iceland",
                        party = "Bright Future",
                        `vote %` = "8.25 (2013)",
                        `Internet daily %` = test_res$estimate[1],
                        `Diff.` = test_res$estimate[1] - test_res$estimate[2],
                        `p` = test_res$p.value,
                        `test` = "Student's t-Test",
                        `source` = "ICENES 2013"))

test_res <-
  t.test(dat_2013$internet_daily[dat_2013$prtvote13 == iceland_parties['Pirate Party']],
       dat_2013$internet_daily[!dat_2013$new_party])

res <- rbind(res, 
             data.frame(country = "Iceland",
                        party = "Pirate Party",
                        `vote %` = "5.1 (2013)",
                        `Internet daily %` = test_res$estimate[1],
                        `Diff.` = test_res$estimate[1] - test_res$estimate[2],
                        `p` = test_res$p.value,
                        `test` = "Student's t-Test",
                        `source` = "ICENES 2013"))

# France

library(haven)
dat_2017 <- 
  read_sav("secondary_data/misc_european_elections/cdsp_FES2017_v2/FES2017.sav")

table(dat_2017$V7)

dat_2017$new_party <- 
  dat_2017$V7 %in% c(`La France insoumise` = 3,
                     `En Marche` = 6,
                     `Debout la France` = 9)

dat_2017$internet_daily <- 
  dat_2017$PA6 %in%  0:1
dat_2017$internet_daily[dat_2017$PA6 > 2] <- NA

test_res <-
  wtd.t.test(dat_2017$internet_daily[dat_2017$new_party],
             dat_2017$internet_daily[!dat_2017$new_party],
             weight = dat_2017$w2[dat_2017$new_party],
             weighty = dat_2017$w2[!dat_2017$new_party])

res <- rbind(res, 
             data.frame(country = "France",
                        party = c("All new parties"),
                        `vote %` = "40.4 (2017)",
                        `Internet daily %` = test_res$additional[2],
                        `Diff.` = test_res$additional[1],
                        `p` = test_res$coefficients[3],
                        `test` = test_res$test,
                        `source` = "FES 2017"))

test_res <-
  wtd.t.test(dat_2017$internet_daily[dat_2017$V7 %in% 3],
             dat_2017$internet_daily[!dat_2017$new_party],
             weight = dat_2017$w2[dat_2017$V7 %in% 3],
             weighty = dat_2017$w2[!dat_2017$new_party])

res <- rbind(res, 
             data.frame(country = "France",
                        party = c("La France insoumise"),
                        `vote %` = "11 (2017)",
                        `Internet daily %` = test_res$additional[2],
                        `Diff.` = test_res$additional[1],
                        `p` = test_res$coefficients[3],
                        `test` = test_res$test,
                        `source` = "FES 2017"))

test_res <-
  wtd.t.test(dat_2017$internet_daily[dat_2017$V7 %in% 6],
             dat_2017$internet_daily[!dat_2017$new_party],
             weight = dat_2017$w2[dat_2017$V7 %in% 6],
             weighty = dat_2017$w2[!dat_2017$new_party])

res <- rbind(res, 
             data.frame(country = "France",
                        party = c("En Marche"),
                        `vote %` = "28.2 (2017)",
                        `Internet daily %` = test_res$additional[2],
                        `Diff.` = test_res$additional[1],
                        `p` = test_res$coefficients[3],
                        `test` = test_res$test,
                        `source` = "FES 2017"))

res <- res[order(as.character(res$country), 
                 as.character(res$`Internet.daily..`)),]
colnames(res)
colnames(res)[3:7] <- c("vote %", "Internet use %", "delta", "p", "test")
res$`Internet use %` <- res$`Internet use %`*100
res$delta <- res$delta *100
res$test <- NULL

res$perc_diff <- 
  res$`Internet use %` /  
  (res$`Internet use %` - res$delta)
  
res$perc_diff <- (res$perc_diff-1)  * 100

# Regression analysis
dat_2014$age <- as.numeric(dat_2014$d11r2)
dat_2014$age[dat_2014$age == 7] <- NA

dat_2014$edu <- as.numeric(dat_2014$d8)
dat_2014$edu[!dat_2014$edu %in% 1:3] <- NA

dat_2014$m5s <- dat_2014$qp2_ees == parties['Five Star Movement']
dat_2014$podemos <- dat_2014$qp2_ees == parties['Podemos']
dat_2014$citizens <- dat_2014$qp2_ees == parties['Citizens']  
dat_2014$demo_left <- dat_2014$qp2_ees == parties['Democratic Left']

res_m5s <- 
  glm(m5s ~ internet_daily + age + edu, 
      family = 'binomial', data = dat_2014[is_italy,])
res_podemos <- 
  glm(podemos ~ internet_daily + age + edu, 
      family = 'binomial', data = dat_2014[is_spain,])
res_citizens <- 
  glm(citizens ~ internet_daily + age + edu, 
      family = 'binomial', data = dat_2014[is_spain,])
res_demo_left <- 
  glm(demo_left ~ internet_daily + age + edu, 
      family = 'binomial', data = dat_2014[is_greece,])
```

### Table 1 (Supplemental online material)

```{r, results='asis'}
require(xtable)
print(xtable(
  res[c("country", "party", "vote %", "Internet use %", "perc_diff",
                   "p", "source")], digits = 1), 
      include.rownames=FALSE, type = "latex")
```


